Tutorials › Multivariate Statistics (Graduate) › HOSVD and HOOI for Tucker Approximation

Tensor and Multiway Models · Tutorial 934 of 975

HOSVD and HOOI for Tucker Approximation

Learn how mode-wise singular vectors produce a Tucker approximation, what its error bounds guarantee, and how alternating subspace updates can improve its fit.

Advanced 12 min read

What You'll Learn

  • Define the truncated HOSVD approximation using leading left singular vectors of tensor unfoldings.
  • Derive the optimal core for fixed factor matrices and the HOSVD residual bound.
  • Compare the HOSVD approximation error with the best approximation at the same multilinear rank.
  • Derive each HOOI update as a leading left singular subspace of a partially projected tensor.
  • Interpret monotone objective improvement without mistaking it for a guarantee of a global optimum.

From Tucker Structure to a Computable Approximation

A Tucker model represents a tensor through a smaller core and one factor matrix for each mode. The central computational question is how to choose those factor subspaces when the desired multilinear rank is specified. The higher-order singular value decomposition (HOSVD) gives a direct construction from the mode unfoldings. Higher-order orthogonal iteration (HOOI) then refines the subspaces by alternating updates.

Let \(\mathcal{X}\in\mathbb{R}^{I_1\times\cdots\times I_N}\) be an order-\(N\) tensor. For mode \(n\), let \(X_{(n)}\) denote its mode-\(n\) unfolding, with mode \(n\) indexing rows and the other modes indexing columns. Unfolding preserves the Frobenius norm. The mode product \(\mathcal{X}\times_n U_n\) applies a matrix \(U_n\) along mode \(n\). Choose target ranks \(R_n\leq I_n\), and let \(U_n\in\mathbb{R}^{I_n\times R_n}\) have orthonormal columns.

Definition: A tensor has multilinear rank at most \((R_1,\ldots,R_N)\) if each unfolding \(X_{(n)}\) has matrix rank at most \(R_n\). Given orthonormal factor matrices \(U_1,\ldots,U_N\), the least-squares Tucker core is \( \mathcal{G}=\mathcal{X}\times_1 U_1^T\cdots\times_N U_N^T \). The associated approximation is \( \widehat{\mathcal{X}}=\mathcal{G}\times_1 U_1\cdots\times_N U_N \).

The factor matrices specify the subspaces; the core stores the coordinates of the tensor in those subspaces. The target ranks are treated as given here. The focus is how to construct and refine the approximation at those ranks.

The Truncated HOSVD

For each mode, compute a singular value decomposition of \(X_{(n)}\). Write its singular values in nonincreasing order as \(\sigma_1^{(n)}\geq\sigma_2^{(n)}\geq\cdots\). Take \(U_n\) to contain \(R_n\) leading left singular vectors. The resulting core and reconstruction are the truncated HOSVD. Define the discarded mode-\(n\) singular-value energy by \( \tau_n^2=\sum_{j>R_n}(\sigma_j^{(n)})^2 \).

Theorem: The truncated HOSVD residual satisfies \( \|\mathcal{X}-\widehat{\mathcal{X}}_{\mathrm{HOSVD}}\|_F^2\leq\sum_{n=1}^N\tau_n^2 \). If \(E_{\mathrm{opt}}\) is the infimum of the Frobenius error over tensors of multilinear rank at most \((R_1,\ldots,R_N)\), then \( E_{\mathrm{opt}}\geq\max_n\tau_n \). Consequently, \( \|\mathcal{X}-\widehat{\mathcal{X}}_{\mathrm{HOSVD}}\|_F\leq\sqrt{N}\,E_{\mathrm{opt}} \).

Proof. Let \(P_n=U_nU_n^T\), the orthogonal projection onto the selected mode-\(n\) subspace. The HOSVD reconstruction is the result of applying all these projections to \(\mathcal{X}\). Projections acting on different modes commute. A telescoping decomposition of the residual is

$$ \mathcal{X}-\mathcal{X}\times_1P_1\cdots\times_NP_N = \sum_{n=1}^N \mathcal{X}\times_1P_1\cdots\times_{n-1}P_{n-1}\times_n(I-P_n). $$

The summands are pairwise orthogonal. Indeed, for two indices \(i<j\), the \(i\)th summand is in the orthogonal complement of the range of \(P_i\), whereas the \(j\)th summand has already been projected into the range of \(P_i\). The squared residual norm is therefore the sum of the squared norms of these summands. Since orthogonal projections in the preceding modes cannot increase the Frobenius norm, the \(n\)th summand has squared norm at most \(\|\mathcal{X}\times_n(I-P_n)\|_F^2\). By the singular value decomposition of \(X_{(n)}\), this last squared norm is exactly \(\tau_n^2\). Summing proves the HOSVD upper bound.

For the lower bound, take any tensor \(\mathcal{Y}\) of multilinear rank at most \((R_1,\ldots,R_N)\). Its mode-\(n\) unfolding has matrix rank at most \(R_n\). The Eckart–Young theorem says that the Frobenius distance from \(X_{(n)}\) to any matrix of rank at most \(R_n\) is at least \(\tau_n\). Since unfolding preserves Frobenius norm, \( \|\mathcal{X}-\mathcal{Y}\|_F\geq\tau_n \) for every \(n\). Taking the infimum over \(\mathcal{Y}\) gives \(E_{\mathrm{opt}}\geq\max_n\tau_n\). Finally, \( \sum_n\tau_n^2\leq N(\max_n\tau_n)^2\leq N E_{\mathrm{opt}}^2 \), which proves the stated comparison. \(\square\)

Worked Example: A Small HOSVD With a Direct Error Check

Consider the order-three \(2\times2\times2\) tensor whose only nonzero entries are \(X_{111}=3\) and \(X_{222}=1\). Use target ranks \((1,1,1)\).

In each unfolding, the two nonzero rows are orthogonal and have norms \(3\) and \(1\). Thus the singular values are \(3,1\), and a leading left singular vector in every mode is \(e_1=(1,0)^T\). The core has the single entry \(G_{111}=X_{111}=3\), so the HOSVD approximation retains only \(3e_1\otimes e_1\otimes e_1\). The residual is \(e_2\otimes e_2\otimes e_2\), with Frobenius norm \(1\).

For each mode, the discarded singular-value energy is \(\tau_n^2=1\). The theorem bounds the residual by \(\sqrt{1+1+1}=\sqrt{3}\), which is valid but not sharp. The unfolding lower bound gives \(E_{\mathrm{opt}}\geq1\); since this HOSVD approximation has error \(1\), it is optimal for these target ranks.

Worked Example: Using the Unfolding Bound to Verify the Optimum

Let the \(2\times2\times2\) tensor have \(X_{111}=2\), \(X_{222}=2\), and all other entries zero. Again use target ranks \((1,1,1)\). In each mode unfolding the two nonzero rows are orthogonal with norm \(2\), so the singular values are \(2,2\). Choose \(e_1\) as the leading singular vector in each mode. The HOSVD approximation is \(2e_1\otimes e_1\otimes e_1\), and its residual is \(2e_2\otimes e_2\otimes e_2\), of norm \(2\).

To establish optimality, note that each discarded tail has \(\tau_n^2=4\). The theorem's unfolding lower bound therefore gives \(E_{\mathrm{opt}}\geq\tau_n=2\). The displayed approximation has error \(2\), so \(E_{\mathrm{opt}}=2\). This lower-bound argument is essential: the observation that the two components cannot both be represented exactly would not, by itself, rule out an approximation with error less than \(2\).

Here the HOSVD upper bound is \(\sqrt{12}\), while the optimal error is \(2\). The general guarantee is a worst-case comparison, not a claim that the HOSVD error always equals that bound.

HOOI: Updating One Subspace at a Time

HOOI seeks to improve the core-subspace approximation by optimizing the factor matrices in turn. The key calculation is the optimal update for one mode when all other factor matrices are held fixed. First, for any fixed orthonormal factors, the best core is the projection of \(\mathcal{X}\) into their coordinates. The residual is orthogonal to the reconstructed tensor, so

$$ \|\mathcal{X}-\widehat{\mathcal{X}}\|_F^2 = \|\mathcal{X}\|_F^2-\|\mathcal{G}\|_F^2. $$

Thus minimizing residual is equivalent to maximizing the squared norm of the projected core. Now fix every factor except \(U_n\), and form the partially projected tensor \( \mathcal{Z}_n=\mathcal{X}\times_{m\ne n}U_m^T \). Its mode-\(n\) unfolding \(Z_{n,(n)}\) has \(I_n\) rows. For a candidate \(U_n\), the resulting core is obtained by multiplying this unfolding by \(U_n^T\) (up to the chosen ordering of columns). Therefore the update solves \( \max_{U_n^TU_n=I}\|U_n^TZ_{n,(n)}\|_F^2 \).

Theorem: An optimal HOOI update for mode \(n\), with all other factors fixed, takes the columns of \(U_n\) to be \(R_n\) leading left singular vectors of \(Z_{n,(n)}\). If singular values tie at the truncation boundary, any maximizing leading subspace is a valid update.

Derivation. Let \(A=Z_{n,(n)}Z_{n,(n)}^T\), a positive semidefinite matrix with eigenvalues \(\lambda_1\geq\lambda_2\geq\cdots\geq0\). Then \( \|U_n^TZ_{n,(n)}\|_F^2=\operatorname{tr}(U_n^TAU_n) \). Write the columns of \(U_n\) in an orthonormal eigenbasis of \(A\). If \(w_j\) is the sum of the squared coefficients assigned to the \(j\)th eigenvector, then \(0\leq w_j\leq1\), \(\sum_jw_j=R_n\), and the trace is \(\sum_j\lambda_jw_j\). This sum is maximized by assigning weight one to the \(R_n\) largest eigenvalues and zero to the rest: moving any positive weight from a smaller eigenvalue to an unfilled larger one cannot decrease the sum. Hence the maximizing columns span a leading eigenspace of \(A\), equivalently a leading left singular subspace of \(Z_{n,(n)}\). \(\square\)

A full HOOI sweep applies this update successively for modes \(1,\ldots,N\), using the most recently updated factors. Each update maximizes the projected-core norm for that mode while keeping the others fixed. It follows that the projected-core norm cannot decrease from one update to the next. The residual norm, by the orthogonal decomposition above, cannot increase. The residual is nonnegative, so its values across updates converge. This monotonicity concerns the objective values; it does not imply that the factors converge or that the final approximation is globally optimal.

Worked Example: A HOOI Update Improves a Rank-One Fit

Consider a \(2\times2\times2\) tensor with \(X_{111}=3\), \(X_{211}=4\), and all other entries zero. Use target ranks \((1,1,1)\). To illustrate one coordinate update, hold \(u_2=u_3=e_1\) and start with \(u_1=e_2\).

The current scalar core is \(u_1^T(X_{:11})=e_2^T(3,4)^T=4\). The current approximation has norm \(4\), and its squared residual is \(3^2=9\), since the tensor has squared norm \(3^2+4^2=25\).

For the mode-1 update, the partially projected mode-1 unfolding reduces to the single column \((3,4)^T\). Its leading left singular vector is \(u_1=(3/5,4/5)^T\), because \(\sqrt{3^2+4^2}=5\). The updated core is \(u_1^T(3,4)^T=5\), so the updated approximation has squared norm \(25\) and residual zero. Indeed, \( \mathcal{X}=(3,4)^T\otimes e_1\otimes e_1 \), so this update exactly recovers the tensor. The example shows the coordinate optimization explicitly; a different starting point, such as the HOSVD factors, may already be optimal.

What the Algorithms Guarantee—and What They Do Not

The HOSVD is a direct construction: one singular value decomposition per unfolding determines all factor subspaces. Its error guarantee compares the residual to discarded energies from each mode and gives a factor-\(\sqrt{N}\) bound relative to the best approximation at the same multilinear rank. HOOI uses repeated singular-value decompositions of partially projected unfoldings to increase the captured core energy. It is common to initialize HOOI with the HOSVD factors, then stop when the residual changes by less than a chosen numerical tolerance.

The HOOI update is an exact solution to a single-factor optimization problem, not to the joint optimization over all factors at once. Since the factors interact through the core, successive coordinate-wise optima need not be a global optimum. Initialization can matter, and an objective that has stopped changing should not be interpreted as proof that no better set of subspaces exists.

Common pitfall: It is incorrect to claim that a rank-constrained approximation must have error at least a particular value merely because it cannot represent two components exactly. Approximation can distribute error across components. To establish a lower bound, use a valid bound such as the Eckart–Young result applied to an unfolding; then compare that bound with the error of a concrete candidate. Likewise, HOOI's monotone decrease in residual does not prove global optimality.

Key takeaway: HOSVD selects leading subspaces independently from the original unfoldings and has a proved error bound. HOOI recomputes each subspace from a tensor projected along the other current factors, guaranteeing nonincreasing residual at each update but not a global solution.

Check Your Understanding

Answer each question using the definitions and guarantees above.

  1. For a fixed set of orthonormal factor matrices, what choice of core minimizes the Frobenius reconstruction error?
  2. Why does the mode-\(n\) discarded singular-value tail give a lower bound on the error of every tensor with multilinear rank at most \((R_1,\ldots,R_N)\)?
  3. In the HOSVD residual proof, why are the terms in the telescoping decomposition pairwise orthogonal?
  4. What matrix is decomposed to obtain a HOOI update for mode \(n\), and which singular vectors are retained?
  5. What does HOOI's monotone residual guarantee establish, and what does it not establish?
Clinical Biostats