PaperlayerView exact revision on arXiv

Source-linked AI summary

Quasi-Newton methods on Grassmannians and multilinear approximations of tensors

Berkant Savas, Lek-Heng Lim

arXiv:0907.2214v3math.OCcs.NAmath.NA

TL;DR

Existing solvers for symmetric multilinear tensor approximation are described as unsatisfactory because they rely on the distinct general-tensor problem. The paper develops Grassmannian quasi-Newton methods and reports speed and robustness improvements over Grassmann-Newton; its BFGS optimality result is coordinate-dependent.

  • Problem

    Existing solvers for symmetric multilinear approximation are described as unsatisfactory because they rely on the distinct general-tensor problem.

    Source: §1.2 · p. 3

  • Method

    The paper develops Grassmannian quasi-Newton methods, using manifold operations represented through matrix operations.

    Source: §5.1 · p. 12, §7.1 · p. 24, §13 · p. 39

  • Results

    The authors report substantial speed and robustness improvements over Grassmann-Newton, while BFGS and L-BFGS converge faster than HOOI in the reported experiments.

    Source: §11.2 · p. 36, §13 · p. 39

  • Takeaways & Limitations

    The paper presents specialized algorithms for symmetric tensor approximation and makes its codes freely available for readers to verify their speed, accuracy, and robustness.

    Source: §13 · p. 39

    The BFGS least-change optimality property holds for local coordinates, not for the global coordinate representation or abstract Hessian operators.

    Source: §6.5 · p. 23

Abstract

from arXiv · show

In this paper we proposed quasi-Newton and limited memory quasi-Newton methods for objective functions defined on Grassmannians or a product of Grassmannians. Specifically we defined BFGS and L-BFGS updates in local and global coordinates on Grassmannians or a product of these. We proved that, when local coordinates are used, our BFGS updates on Grassmannians share the same optimality property as the usual BFGS updates on Euclidean spaces. When applied to the best multilinear rank approximation problem for general and symmetric tensors, our approach yields fast, robust, and accurate algorithms that exploit the special Grassmannian structure of the respective problems, and which work on tensors of large dimensions and arbitrarily high order. Extensive numerical experiments are included to substantiate our claims.

1. Introduction.

The paper develops BFGS and L-BFGS methods that respect Grassmannian geometry, including a local-coordinate update with a Euclidean-style optimality property. It applies these methods to general and symmetric tensor approximation, reducing the problems to Grassmannian optimization and reporting improvements over existing approaches.

1.1. Quasi-Newton and limited memory quasi-Newton algorithms on Grassmannians.

  • The algorithms define BFGS and L-BFGS on Grassmannians and products of Grassmannians using standard numerical linear algebra operations.Geodesic movement and parallel transport are represented by matrix operations rather than numerical solutions of defining differential equations.Source: §1.1 · p. 1
  • Local-coordinate BFGS retains the Euclidean update’s optimality property, while the corresponding property does not hold for global-coordinate BFGS.The local update is the best possible Hessian-approximation update satisfying the secant equations and preserving symmetric positive definiteness.Source: §1.1 · p. 2
  • For products of Grassmannians, local-coordinate BFGS has unit cost O(kr(n −r)), compared with O(krn) for global-coordinate BFGS.For a single Grassmannian, the corresponding costs are r(n−r) and nr.Source: §1.1 · pp. 1–2

1.2. Multilinear approximations of tensors and symmetric tensors.

  • The tensor problems are reduced to maximization over products of Grassmannians in the general case and a single Grassmannian in the symmetric case.The reduction uses invariance under right orthogonal transformations; the general formulation eliminates the core tensor from the objective.Source: §1.2 · p. 3
  • The paper addresses symmetric tensor approximation directly, rather than imposing symmetry after solving a general tensor problem.The authors note that using the general framework without accounting for symmetry incurs a k-fold computational cost relative to the symmetric formulation.Source: §1.2 · p. 3
  • The authors report guaranteed convergence to stationary points and better speed and robustness than Grassmann-Newton, with accuracy comparable to it and better than Gauss-Seidel strategies in many cases.The paper gives examples including order-10 tensor approximation and reports that large-problem L-BFGS can outperform HOOI.Source: §1.2 · p. 4

2. General and symmetric tensors.

The section introduces general and symmetric tensors and illustrates why symmetric tensors arise in applications. It establishes notation for tensor and symmetric-tensor representations.

2. General and symmetric tensors.

  • A tensor of order k is represented as a k-dimensional array, while a symmetric tensor has entries invariant under permutations of its indices.The section assumes fixed bases and sometimes uses the order-3 case for notational simplicity.Source: §2 · p. 5
  • Higher derivatives of smooth functions and moments and cumulants of random vectors provide examples of symmetric tensors.For a scalar random variable, the first four cumulants are expectation, variance, skewness, and kurtosis.Source: §2 · pp. 5–6

3. Multilinear transformation and multilinear rank.

The section defines multilinear tensor transformations and their associated multilinear ranks, then formulates best-approximation problems as Grassmannian optimization. This reduction supports dimension reduction for both general and symmetric tensors.

3. Multilinear transformation and multilinear rank.

  • Multilinear rank records the ranks of a tensor’s mode-wise matrix unfoldings and is generally a tuple rather than one number.A rank-one tensor has multilinear rank (1, 1, 1).Source: §3 · p. 7
  • Symmetric multilinear rank is the common mode rank for a symmetric tensor, and this paper uses that meaning whenever it refers to symmetric tensor rank.The paper distinguishes it from rank defined by the minimum number of rank-one terms; it excludes that other approximation problem because it may be ill-posed and may have no solution.Source: §3 · p. 8

3.1. Multilinear approximation as maximization over Grassmannians.

  • For general tensors, orthonormal factor matrices yield a best multilinear-rank approximation objective on a product of Grassmannians, with the core tensor eliminated.The objective is invariant to right multiplication of each factor by an orthogonal matrix.Source: §3.1 · p. 9
  • For symmetric tensors, the corresponding best symmetric multilinear-rank approximation is formulated as optimization on one Grassmannian.The objective is invariant to right orthogonal transformations of its factor matrix.Source: §3.1 · p. 10
  • A tensor approximation reduces storage from O(n^k) to O(knr + rk), which can yield orders-of-magnitude size reduction when r is sufficiently smaller than n.The passage gives r = O(n^(1/k)) as an example of a rank scale where this may occur.Source: §3.1 · p. 10

4. Optimization in Euclidean space and on Riemannian manifolds.

The section reviews Euclidean BFGS before explaining which operations must change for optimization on a Riemannian manifold. On manifolds, gradients and Hessians act on tangent spaces, iterates move along geodesics, and vectors and operators require transport between points.

4.1. BFGS updates in Euclidean Space.

  • Euclidean BFGS updates a Hessian approximation from one iterate to the next instead of recomputing the exact Hessian at every step.The search direction is obtained from the current iterate and gradient, and the step length is typically chosen by a line search satisfying Wolfe or Goldstein conditions.Source: §4.1 · p. 11

4.2. Quasi-Newton methods on a Riemannian manifold.

  • On a manifold, the gradient lies in the tangent space at the current point, and differentiating a global-coordinate expression does not generally give the correct gradient or Hessian.The Hessian is a linear transformation of that point’s tangent space to itself.Source: §4.2 · pp. 11–12
  • A Euclidean step may leave the manifold, so updates move along geodesics; parallel transport aligns tangent vectors and operators defined at different points.The update’s terms must be defined at the same point for the BFGS operations to be well-defined.Source: §4.2 · p. 12

5. Grassmann geodesics and parallel transport of vectors.

Grassmannian geodesics and parallel transport can be expressed with matrix operations, enabling algorithms to use standard matrix arithmetic. The section distinguishes global and local tangent coordinates and explains how local coordinates simplify transport at the cost of maintaining a tangent-space basis.

5. Grassmann geodesics and parallel transport of vectors.

  • The section frames geodesics and parallel transport as matrix operations, including for products of Grassmannians.It identifies following a geodesic and transporting a tangent along it as the two operations needed here.Source: §5.1 · p. 12

5.2. Grassmannians in terms of matrices.

  • A point is represented by an orthonormal matrix from its equivalence class, while a tangent Δ can be written as X⊥D in global and local coordinates, respectively.The local representation D has r(n−r) degrees of freedom.Source: §5.2 · p. 13

5.3. Geodesics.

  • A geodesic from X in tangent direction Δ is represented using the thin SVD of Δ; retaining its final V factor matters for representation consistency in quasi-Newton Grassmann algorithms.Omitting that factor traces the same manifold path with a different matrix representation.Source: §5.3 · pp. 13–14

5.4. Parallel transport in global and local coordinates.

  • Local coordinates keep a tangent’s matrix representation constant along a geodesic, avoiding explicit transport of tangents and Hessians while requiring transport of the basis X⊥.For the stated operations, tangent and Hessian transport costs are O(n^2r) and O(n^3r^2); updating X⊥ costs at most O(n^3) initially and O(n^2(n−r)) per later iteration.Source: §5.4 · pp. 14–15

6. Quasi-Newton methods with BFGS updates on a Grassmannian.

This section develops BFGS updates for Grassmannians and products of Grassmannians in global and local coordinates, including tensor representations, and shows that local tangent-space coordinates recover the Euclidean BFGS update and its least-change optimality property. It also introduces Euclidean L-BFGS background and notes alternative manifold optimization operations whose convergence effects are unclear.

6.1. Computations in global coordinates.

  • Global-coordinate updates form gradient, step, and gradient-difference quantities in tangent spaces, transport the Hessian between iterates, and then apply BFGS.The transported Hessian is represented by multiplying by transport matrices on both sides before the update.Source: §6.1 · pp. 16–17

6.2. Computations in local coordinates.

  • In local coordinates, transporting the orthogonal-complement basis keeps tangent coordinates constant, so the Hessian representation stays fixed and BFGS uses the Euclidean update.The update avoids explicit Hessian parallel transport, unlike the global-coordinate update.Source: §6.2 · p. 17

6.3. BFGS update in tensor form.

  • The tensor formulation keeps tangent steps and gradient differences as matrices and expresses Hessian actions and BFGS updates through tensor contractions and products.In local coordinates, the tensor update is simpler because the Hessian does not need parallel transport.Source: §6.3 · p. 18

6.4. BFGS update on a product of Grassmannians.

  • For a product of Grassmannians, the Hessian is a block tensor operator whose blocks map between tangent spaces, and BFGS extends by applying that operator and using block tensor products.The blocks can have different, nonconforming dimensions.Source: §6.4 · pp. 18–19

6.5. Optimality of BFGS on Grassmannians.

  • In local tangent coordinates, the BFGS update is the Euclidean rank-2 update and is optimal when expressed in the transported tangent basis.The local Hessian representation remains unchanged along a geodesic when its tangent basis is transported.Source: §6.5 · pp. 19–21
  • Theorem 6.6 establishes that local-coordinate BFGS on a product of Grassmannians is the least-change Hessian update satisfying the secant equations.The result follows because the local Hessian representation is constant along geodesics when tangent-space bases are tracked.Source: §6.5 · p. 23
  • The optimality result is coordinate-dependent: least-change optimality does not hold for global-coordinate representations or abstract Hessian operators.Source: §6.5 · p. 23
  • Retractions and vector transport are alternatives to geodesic movement and parallel transport, but the convergence effects of general vector transport are unclear.Source: §6.6 · p. 23

6.6. Other alternatives.

  • The section introduces Euclidean L-BFGS updates as background for the Grassmann variant, citing and [52, Chapter 7].It states that the Euclidean BFGS update has a compact representation, but this passage does not provide the formula.Source: §7 · p. 23

7. Limited memory BFGS.

The section adapts limited-memory BFGS to Grassmannians by retaining a small history of update vectors and transporting those vectors between tangent spaces. It also describes implementation trade-offs for large, sparse problems.

7. Limited memory BFGS.

  • Keeping only the m most recent update pairs and using a scaled identity for the initial Hessian makes each iteration’s storage and computation a small fraction of regular BFGS.The text reports satisfactory results often use 5 ≤ m ≤ 20, including for large problems.Source: §7 · p. 24

7.1. Limited memory BFGS on Grassmannians.

  • Grassmann L-BFGS transports stored tangent vectors between iterations, while inner products between transported vectors remain unchanged and the tangent space is invariant under the update operator.The transported operator remains positive definite and maps tangent vectors to tangent vectors.Source: §7.1 · pp. 24–25
  • For large, sparse problems, explicitly transporting the 2(m−1) stored vectors can be more economical than storing and updating the local-coordinate basis X⊥.The text notes that m is usually small compared with Hessian dimensions.Source: §7.1 · p. 25

8. Quasi-Newton methods for the best multilinear rank approximation

The section formulates best multilinear rank approximation as optimization over Grassmannians and gives gradient and Hessian expressions for these objectives. It extends the expressions to order-k tensors using local tangent coordinates.

8. Quasi-Newton methods for the best multilinear rank approximation

  • For a three-factor tensor, the objective maximizes the squared Frobenius norm of the contracted core subject to orthonormality of the factor matrices.The factors represent points on a product of three Grassmannians, and the Grassmann gradient has one component for each factor.Source: §8 · p. 25
  • The global-coordinate gradient components are obtained by projecting each Euclidean derivative onto its factor’s tangent space.The projections use ΠX = I − XXT, with corresponding expressions for the other factors.Source: §8 · p. 25

8.1. General expression for Grassmann gradients and Hessians.

  • For order-k tensors, the local-coordinate gradient is a k-tuple whose ith component is the contraction ⟨Bi, F⟩−i.Each component is an (ni−ri) × ri matrix representing a tangent at Xi.Source: §8.1 · p. 26
  • The order-k Grassmann Hessian is given in implementable form, with diagonal blocks as Sylvester operators and separate off-diagonal block operators.The text states the result directly and refers readers to for its derivation.Source: §8.1 · p. 27

9. Best multilinear rank approximation of a symmetric tensor.

For noisy symmetric tensors, the section formulates low-rank approximation on one Grassmannian and derives the gradient and Hessian, including for arbitrary tensor order.

9. Best multilinear rank approximation of a symmetric tensor.

  • Symmetry forces all mode ranks to equal r, and HOSVD yields S=(X,...,X)·C; the noisy approximation is optimized on one Grassmannian.This uses less storage and computation than the general problem of the same dimensions.Source: §9 · pp. 27–28

9.1. The symmetric Grassmann gradient.

  • The symmetric Grassmann gradient reuses the general-case tensor-product terms; symmetry makes the modewise expressions equal, so one expression suffices.The derivation uses the first equivalent modewise expression without loss of generality.Source: §9 · p. 28, §9.1 · p. 28, §9.3 · p. 29

9.2. The symmetric Grassmann Hessian.

  • For symmetric 3-tensors, 15 Taylor second-order terms reduce by symmetry to three Hessian terms, including a Grassmann-specific term.The Grassmann-specific term arises because the objective is defined on a Grassmannian.Source: §9.1 · p. 28, §9.2 · p. 29

9.3. General expression for Grassmann gradients and Hessians for a symmetric tensor.

  • For symmetric k-tensors, symmetry makes the diagonal blocks identical and the off-diagonal blocks identical; the complete Hessian combines both.The diagonal and off-diagonal block operators are summed as H = Hdiag + Hoff-diag.Source: §9.3 · pp. 29–30

9.4. Matricizing the Hessian operator.

  • The Hessian operator is matricized into H1+H2+H3; local coordinates replace ΠX with X⊥ before applying the same matricization procedure.Tensor matricization maps the contractions into matrix-vector products.Source: §9.4 · p. 30

10. Examples.

The examples demonstrate gradient calculations and geodesic steps on Grassmannians, including a worked objective increase for a 3×3×3 tensor.

Examples

  • The example computes gradients globally and locally; local coordinates parallel-transport chosen bases, so their representation depends on the basis choice.The orthogonal-complement bases are parallel transported along the geodesic.Source: §10 · pp. 31–32
  • 40.5→≈45.5625 after a π/4 step, the expected objective increase in the worked 3×3×3 tensor example.These are the objective values at the starting point and the new point, respectively.Source: §10 · p. 33
  • For Gr(3,1)=S², the procedure chooses an ascent tangent, moves along a great-circle geodesic, then repeats at the new point.For rank-1 approximation of a 3×3×3 tensor, optimization takes place on a product of three spheres.Source: §10 · p. 33

11. Numerical experiments and computational complexity.

Experiments compare Grassmann BFGS and L-BFGS with Newton-Grassmann and HOOI across general, symmetric, higher-order, and mixed-order tensor problems. The section also analyzes coordinate costs and the memory boundary imposed by explicitly stored tensors.

11.1. Initialization and stopping condition.

  • Experiments use a relative-gradient stopping ratio typically near 10^-13, with HOOI iterations used to initialize points and make the initial Hessian negative definite.BFGS starts from an exact Hessian or scaled identity; L-BFGS uses a scaled identity, with m between 5 and 30.Source: §11 · p. 33, §11.1 · p. 34

11.2. Experimental results.

  • BFGS converges superlinearly on a 20 × 20 × 20 general tensor, while L-BFGS with m=20 handles a 200 × 200 × 200 case; varying m from 10 to 30 yields indistinguishable runs.For the 50 × 50 × 50 rank-(20, 20, 20) case, m=5 generally performs slightly worse than m=10, 15, 20, 25, or 30.Source: §11.2 · p. 34, Figure 11.1 · p. 34, Figure 11.2 · p. 35
  • For symmetric 50 × 50 × 50 approximation, exact-Hessian BFGS converges much faster than scaled-identity BFGS, nearly matching Newton-Grassmann.The figure also reports a 100 × 100 × 100 symmetric tensor approximated at rank 20, with m=10.Source: §11.2 · p. 35, Figure 11.3 · p. 35
  • BFGS converges rapidly near local minima on order-10 and mixed-order objectives; experiments also include general and symmetric 4-tensor approximations.The mixed objective combines order-2, order-3, and order-4 tensors; the examples report fast convergence near a local minimizer.Source: §11.2 · pp. 35–36, Figure 11.5 · p. 36, Figure 11.4 · p. 36
  • The convergence plots show BFGS and L-BFGS converge faster than HOOI, and MATLAB implementations give them the shortest runtime; among currently available methods, these quasi-Newton methods have best computational efficiency.Source: §11.2 · p. 36, §11.3 · p. 38

11.3. Computational complexity, curse of dimensionality, and conver-

  • For general 3-tensors, local coordinates reduce the Hessian dimension from 3nr × 3nr to 3(n−r)r × 3(n−r)r; the symmetric problem differs in complexity only by a constant factor.Table 11.1 compares computational complexity across global-coordinate BFGS, local-coordinate BFGS, and L-BFGS.Source: §11.3 · pp. 37–38, Table 11.1 · p. 37
  • Explicitly stored order-k tensors require n^k entries, and the methods apply while the given tensor fits in memory; compact or functional representations can supply the needed derivatives directly.The stated additional memory is O(n^2r^2k^2) for BFGS and O(2nrkm) for L-BFGS.Source: §11.3 · p. 38

12. Related work.

Related tensor-approximation methods differ in how they converge and use tensor information; the paper presents Grassmannian quasi-Newton algorithms as another approach.

12. Related work.

  • ALS-type convergence appears to depend on the decay of multilinear singular values, though the exact dependence is unclear.The cited methods include ALS and multigrid-accelerated ALS.Source: §12 · p. 38
  • Cross methods retain only a fraction of the tensor and assume approximations within specified bounds and ranks, unlike this paper’s use of all given information.The paper describes its approach as making no such assumptions.Source: §12 · p. 39
  • The proposed BFGS and L-BFGS methods respect the Riemannian metric structure and use standard matrix operations, with local-coordinate updates extending BFGS optimality to Grassmannians and their products.The paper identifies local coordinates as a novelty in manifold optimization.Source: §13 · p. 39

13. Conclusion.

The paper applies its algorithms to general and symmetric tensor approximation, reports improvements over Grassmann-Newton in speed and robustness, and presents symmetry-specialized algorithms.

13. Conclusion.

  • The algorithms address best multilinear rank approximation for general tensors and an analogous, but very different, problem for symmetric tensors.The conclusion also reports substantial speed and robustness improvements over the Grassmann-Newton algorithm.Source: §13 · p. 39
  • Specialized algorithms account for symmetry in multilinear approximation of symmetric tensors and related problems.Source: §13 · p. 39
  • The authors made their codes freely available so readers can verify the algorithms’ speed, accuracy, and robustness.Source: §13 · p. 39

Appendix A. Notation for tensor contractions.

The appendix introduces contracted tensor-product notation, including subscripts and semicolons to indicate which indices are contracted.

Appendix A. Notation for tensor contractions.

  • For third-order tensors, the appendix introduces notation for contracted products.Source: Appendix A · p. 39
  • Subscripts such as 1 or 1,2 indicate contraction over the first index or the first two indices, respectively.Source: Appendix A · p. 39
  • A semicolon separates the index lists when contraction is performed over different indices, and superscript-style negative-index notation can denote contraction over all but specified indices.For example, ⟨A, B⟩1,2 = ⟨A, B⟩−3 and ⟨A, B⟩1 = ⟨A, B⟩−(2,3).Source: Appendix A · pp. 39–40
Loading 0907.2214v3…