---
abstract: |
    In many applications such as data compression, imaging or genomic data analysis, it is important to approximate a given tensor by a tensor that is sparsely representable. For matrices, i.e. $2$-tensors, such a representation can be obtained via the singular value decomposition, which allows to compute best rank $k$-approximations. For very big matrices a low rank approximation using SVD is not computationally feasible. In this case different approximations are available. It seems that variants of the CUR-decomposition are most suitable.

    For $d$-mode tensors ${\mathcal{T}}\in \otimes_{i=1}^d {\mathbb{R}}^{n_i}$, with $d>2$, many generalizations of the singular value decomposition have been proposed to obtain low tensor rank decompositions. The most appropriate approximation seems to be best $(r_1,\ldots,r_d)$-approximation, which maximizes the $\ell_2$ norm of the projection of ${\mathcal{T}}$ on $\otimes_{i=1}^d {\mathbf{U}}_i$, where ${\mathbf{U}}_i$ is an $r_i$-dimensional subspace ${\mathbb{R}}^{n_i}$. One of the most common methods is the *alternating maximization method* (AMM). It is obtained by maximizing on one subspace ${\mathbf{U}}_i$, while keeping all other fixed, and alternating the procedure repeatedly for $i=1,\ldots,d$. Usually, AMM will converge to a local best approximation. This approximation is a fixed point of a corresponding map on Grassmannians. We suggest a Newton method for finding the corresponding fixed point. We also discuss variants of CUR-approximation method for tensors.

    The first part of the paper is a survey on low rank approximation of tensors. The second new part of this paper is a new Newton method for best $(r_1,\ldots,r_d)$-approximation. We compare numerically different approximation methods.
author:
- |
    Shmuel Friedlandand Venu Tammali\
    Department of Mathematics, Statistics and Computer Science,\
    University of Illinois at Chicago,\
    Chicago, Illinois 60607-7045, USA,\
    E-mail: `friedlan@uic.edu, vtamma2@uic.edu`
date: 'December 10, 2014'
title: 'Low-rank approximation of tensors'
---

[**2000 Mathematics Subject Classification.**]{} 14M15, 15A18, 15A69, 65H10, 65K10.

[**Key words**]{} Tensor, best rank one approximation, best $(r_1,\ldots,r_d)$-approximation, sampling, alternating maximization method, singular value decomposition, Grassmann manifold, fixed point, Newton method.

Introduction {#S:intro}
============

Let ${\mathbb{R}}$ be the field of real numbers. Denote by ${\mathbb{R}}^{\mathbf{n}}={\mathbb{R}}^{n_1\times\ldots\times n_d}:=\otimes_{i=1}^d {\mathbb{R}}^{n_j}$, where ${\mathbf{n}}=(n_1,\ldots,n_d)$, the tensor products of ${\mathbb{R}}^{n_1},\ldots,{\mathbb{R}}^{n_d}$. ${\mathcal{T}}=[t_{i_1,\ldots, i_d}]\in {\mathbb{R}}^{\mathbf{n}}$ is called a *$d$-mode tensor*. Note that the number of coordinates of ${\mathcal{T}}$ is $N=n_1\ldots n_d$. A tensor ${\mathcal{T}}$ is called a *sparsely representable tensor* if it can represented with a number of coordinates that is much smaller than $N$.

Apart from sparse matrices, the best known example of a sparsely representable $2$-tensor is a low rank approximation of a matrix $A\in {\mathbb{R}}^{n_1\times n_2}$. A rank $k$-approximation of $A$ is given by $A_{\textrm{appr}}:=\sum_{i=1}^k {\mathbf{u}}_i{\mathbf{v}}_i{^\top}$, which can be identified with $\sum_{i=1}^k {\mathbf{u}}_i\otimes {\mathbf{v}}_i$. To store $A_{\textrm{appr}}$ we need only the $2k$ vectors ${\mathbf{u}}_1,\ldots,{\mathbf{u}}_k\in{\mathbb{R}}^{n_1},\;{\mathbf{v}}_1,\ldots,{\mathbf{v}}_k\in{\mathbb{R}}^{n_2}$. A best rank $k$-approximation of $A\in {\mathbb{R}}^{n_1\times n_2}$ can be computed via the *singular value decomposition*, abbreviated here as SVD, [@GolV96]. Recall that if $A$ is a real symmetric matrix, then the best rank $k$-approximation must be symmetric, and is determined by the spectral decomposition of $A$.

The computation of the SVD requires $\mathcal {O}(n_1n_2^2)$ operations and at least $\mathcal{O} (n_1 n_2)$ storage, assuming that $n_2\le n_1$. Thus, if the dimensions $n_1$ and $n_2$ are very large, then the computation of the SVD is often infeasible. In this case other type of low rank approximations are considered, see e.g. [@AM01; @DV06; @DKM06; @FriKNZ06; @FMMN11; @FKV04; @GTZ97].

For $d$-tensors with $d>2$ the situation is rather unsatisfactory. It is a major theoretical and computational problem to formulate good generalizations of low rank approximation for tensors and to give efficient algorithms to compute these approximations, see e.g. [@deLdV00; @LMV00; @ES09; @FMMN11; @FMPS13; @Kho; @KB09; @MMD06; @Ose11; @OsT09; @ZhaG01].

We now discuss briefly the main ideas of the approximation methods for tensors discussed in this paper. We need to introduce (mostly) standard notation for tensors. Let $[n]:=\{1,\ldots,n\}$ for $n\in{\mathbb{N}}$. For ${\mathbf{x}}_i:=(x_{1,i},\ldots,x_{n_i,i}){^\top}\in {\mathbb{R}}^{n_i},i\in[d]$, the tensor $\otimes_{i\in[d]}{\mathbf{x}}_i={\mathbf{x}}_1\otimes\cdots\otimes {\mathbf{x}}_d={\mathcal{X}}=[x_{j_1,\ldots,j_d}]\in{\mathbb{R}}^{\mathbf{n}}$ is called a decomposable tensor, or rank one tensor if ${\mathbf{x}}_i\ne{\mathbf{0}}$ for $i\in [d]$. That is, $x_{j_1,\ldots,j_d}=x_{j_1,1}\cdots x_{j_d,d}$ for $j_i\in[n_i], i\in[d]$. Let ${{\langle}{\mathbf{x}}_i,{\mathbf{y}}_i{\rangle}}_i:={\mathbf{y}}_i{^\top}{\mathbf{x}}_i$ be the standard inner product on ${\mathbb{R}}^{n_i}$ for $i\in[d]$. Assume that ${\mathcal{S}}=[s_{j_1,\ldots,j_d}]$ and ${\mathcal{T}}=[t_{j_1,\ldots,j_d}]$ are two given tensors in ${\mathbb{R}}^{\mathbf{n}}$. Then ${{\langle}{\mathcal{S}},{\mathcal{T}}{\rangle}}:=\sum_{j_i\in[n_i],i\in[d]} s_{j_1,\ldots,j_d}t_{j_1,\ldots,j_d}$ is the standard inner product on ${\mathbb{R}}^{\mathbf{n}}$. Note that $$\begin{aligned}
 &&{{\langle}\otimes_{i\in[d]}{\mathbf{x}}_i,\otimes_{i\in[d]}{\mathbf{y}}_i{\rangle}}=\prod_{i\in[d]}{{\langle}{\mathbf{x}}_i,{\mathbf{y}}_i{\rangle}}_i,\\
 &&{{\langle}{\mathcal{T}},\otimes_{i\in[d]}{\mathbf{x}}_i{\rangle}}=\sum_{j_i\in[n_i],i\in[d]} t_{j_1,\ldots,j_d}x_{j_1,1}\cdots x_{j_d,d}.
 \end{aligned}$$ The norm $\|{\mathcal{T}}\|:=\sqrt{{{\langle}{\mathcal{T}},{\mathcal{T}}{\rangle}}}$ is called the Hilbert-Schmidt norm. (For matrices, i.e. $d=2$, it is called the Frobenius norm.)

Let $I=\{1\le i_1<\cdots <i_l\le d\}\subset [d]$. Assume that ${\mathcal{X}}=[x_{j_{i_1},\cdots,j_{i_l}}]\in \otimes_{k\in[l]} {\mathbb{R}}^{n_{i_k}}$. Then the contraction ${\mathcal{T}}\times {\mathcal{X}}$ on the set of indices $I$ is given by: $${\mathcal{T}}\times {\mathcal{X}}=\sum_{j_{i_k}\in [n_{i_k}], k\in[l]} t_{j_1,\ldots,j_d} x_{j_{i_1},\ldots,j_{i_l}}\in\otimes_{p\in[d]\setminus I} {\mathbb{R}}^{n_p}.$$

Assume that ${\mathbf{U}}_i\subset {\mathbb{R}}^{n_i}$ is a subspace of dimension $r_i$ with an orthonormal basis ${\mathbf{u}}_{1,i},\ldots,{\mathbf{u}}_{r_i,i}$ for $i\in[d]$. Let ${\mathbf{U}}:=\otimes_{i=1}^d {\mathbf{U}}_i\subset {\mathbb{R}}^{\mathbf{n}}$. Then $\otimes_{i=1}^d {\mathbf{u}}_{j_i,i}$, where $j_i\in[n_i],i\in[d]$, is an orthonormal basis in ${\mathbf{U}}$. We are approximating ${\mathcal{T}}\in{\mathbb{R}}^{n_1\times \cdots\times n_d}$ by a tensor $$\label{Sapprox}
 {\mathcal{S}}=\sum_{j_i\in[r_i],i\in[d]} s_{j_1,\ldots,j_d} {\mathbf{u}}_{j_1,1}\otimes\cdots  \otimes{\mathbf{u}}_{j_d,d} \in {\mathbb{R}}^{\mathbf{n}}$$ The tensor ${\mathcal{S}}'=[s_{j_1,\ldots,j_d}]\in {\mathbb{R}}^{r_1\times \cdots\times r_d}$ is the *core tensor* corresponding to ${\mathcal{S}}$ in the terminology of [@Tuc66].

There are two major problems: The first one is how to choose the subspaces ${\mathbf{U}}_1,\ldots,{\mathbf{U}}_d$. The second one is the choice of the core tensor ${\mathcal{S}}'$. Suppose we already made the choice of ${\mathbf{U}}_1,\ldots,{\mathbf{U}}_d$. Then ${\mathcal{S}}=P_{{\mathbf{U}}}({\mathcal{T}})$ is the orthogonal projection of ${\mathcal{T}}$ on ${\mathbf{U}}$: $$\label{PUfor}
 P _{\otimes_{i\in[d]}{\mathbf{U}}_i}({\mathcal{T}})=\sum_{j_i\in[r_i],i\in[d]} {{\langle}{\mathcal{T}},\otimes_{i\in[d]} {\mathbf{u}}_{j_i,i}{\rangle}} \otimes_{i\in[d]} {\mathbf{u}}_{j_i,i}.$$ If the dimensions of $n_1,\ldots,n_d$ are not too big, then this projection can be explicitly carried out. If the dimension $n_1,\ldots,n_d$ are too big to compute the above projection, then one needs to introduce other approximations. That is, one needs to compute the core tensor ${\mathcal{S}}'$ appearing in accordingly. The papers [@AM01; @DV06; @DKM06; @FriKNZ06; @FMMN11; @FKV04; @GTZ97; @Kho; @MMD06; @Ose11; @OsT09] essentially choose ${\mathcal{S}}'$ in a particular way.

We now assume that the computation of $P_{{\mathbf{U}}}({\mathcal{T}})$ is feasible. Recall that $$\label{normprojfor}
 \|P_{\otimes_{i\in[d]}{\mathbf{U}}_i}({\mathcal{T}})\|^2=\sum_{j_i\in [r_i],i\in[d]}|{{\langle}{\mathcal{T}},\otimes_{i=1}^d {\mathbf{u}}_{j_i,i}{\rangle}}|^2.$$ The best ${\mathbf{r}}$-approximation of ${\mathcal{T}}$, where ${\mathbf{r}}= (r_1,\ldots,r_d)$, in Hilbert-Schmidt norm is the solution of the minimal problem: $$\label{bestapproxprb}
 \min_{{\mathbf{U}}_i,\dim {\mathbf{U}}_i=r_i, i\in[d]}\min_{{\mathcal{X}}\in \otimes_{i\in[d]}^d{\mathbf{U}}_i} \|{\mathcal{T}}-{\mathcal{X}}\|.$$ This problem is equivalent to the following maximum $$\label{maxnormprb}
 \max_{{\mathbf{U}}_i,\dim {\mathbf{U}}_i=r_i, i\in[d]} \|P_{\otimes_{i\in[d]}{\mathbf{U}}_i}({\mathcal{T}})\|^2.$$

The standard *alternating maximization method*, denoted by AMM, for solving is to solve the maximum problem, where all but the subspace ${\mathbf{U}}_i$ is fixed. Then this maximum problem is equivalent to finding an $r_i$-dimensional subspace of ${\mathbf{U}}_i$ containing the $r_i$ biggest eigenvalues of a corresponding nonnegative definite matrix $A_i({\mathbf{U}}_1,\ldots,{\mathbf{U}}_{i-1},{\mathbf{U}}_{i+1},\ldots,{\mathbf{U}}_d)\in
 {\mathrm{S}}_{n_i}$. Alternating between ${\mathbf{U}}_1,{\mathbf{U}}_2,\ldots,{\mathbf{U}}_d$ we obtain a nondecreasing sequence of norms of projections which converges to $v$. Usually, $v$ is a critical value of $\|P_{\otimes_{i\in[d]}{\mathbf{U}}_i}({\mathcal{T}})\|$. See [@LMV00] for details.

Assume that $r_i=1$ for $i\in [d]$. Then $\dim{\mathbf{U}}_i=1$ for $i\in[d]$. In this case the minimal problem is called a best rank one approximation of ${\mathcal{T}}$. For $d=2$ a best rank one approximation of a matrix ${\mathcal{T}}=T\in{\mathbb{R}}^{n_1\times n_2}$ is accomplished by the first left and right singular vectors and the corresponding maximal singular value $\sigma_1(T)$. The complexity of this computation is ${\mathcal{O}}(n_1n_2)$ [@GolV96]. Recall that the maximum is equal to $\sigma_1(T)$, which is also called the spectral norm $\|T\|_2$. For $d>2$ the maximum is called the spectral norm of ${\mathcal{T}}$, and denoted by $\|{\mathcal{T}}\|_\sigma$. The fundamental result of Hillar-Lim [@HL13] states that the computation of $\|{\mathcal{T}}\|_\sigma$ is NP-hard in general. Hence the computation of best ${\mathbf{r}}$-approximation is NP-hard in general.

Denote by ${\mathop{\mathrm{Gr}}\nolimits}(r,{\mathbb{R}}^n)$ the variety of all $r$-dimensional subspaces in ${\mathbb{R}}^n$, which is called Grassmannian or Grassmann manifold. Let $${\mathbf{1}}_d:=(\underbrace{1,\ldots,1}_d), \quad {\mathop{\mathrm{Gr}}\nolimits}({\mathbf{r}},{\mathbf{n}}):={\mathop{\mathrm{Gr}}\nolimits}(r_1,n_1)\times \cdots \times {\mathop{\mathrm{Gr}}\nolimits}(r_d,n_d).$$ Usually, the AMM for best ${\mathbf{r}}$-approximation of ${\mathcal{T}}$ will converge to a fixed point of a corresponding map ${\mathbf{F}}_{\mathcal{T}}:{\mathop{\mathrm{Gr}}\nolimits}({\mathbf{r}},{\mathbf{n}})\to {\mathop{\mathrm{Gr}}\nolimits}({\mathbf{r}},{\mathbf{n}})$. This observation enables us to give a new Newton method for finding a best ${\mathbf{r}}$-approximation to ${\mathcal{T}}$. For best rank one approximation the map ${\mathbf{F}}_{\mathcal{T}}$ and the corresponding Newton method was stated in [@FMPS11].

This paper consists of two parts. The first part surveys a number of common methods for low rank approximation methods of matrices and tensors. We did not cover all existing methods here. We were concerned mainly with the methods that the first author and his collaborators were studying, and closely related methods. The second part of this paper is a new contribution to Newton algorithms related to best ${\mathbf{r}}$-approximations. These algorithms are different from the ones given in [@ZhaG01; @ES09; @SL10]. Our Newton algorithms are based on finding the fixed points corresponding to the map induced by the AMM. In general its known that for big size problem, where each $n_i$ is big for $i\in[d]$ and $d\ge 3$, Newton methods are not efficient. The computation associate the matrix of derivatives (Jacobian) is too expensive in computation and time. In this case AMM or MAMM (modified AMM) are much more cost effective. This well known fact is demonstrated in our simulations.

We now briefly summarize the contents of this paper. In §\[S:SVD\] we review the well known facts of singular value decomposition (SVD) and its use for best rank $k$-approximation of matrices. For large matrices approximation methods using SVD are computationally unfeasible. Section \[S:lowrankapmat\] discusses a number of approximation methods of matrices which do not use SVD. The common feature of these methods is sampling of rows, or columns, or both to find a low rank approximation. The basic observation in §\[Sub:rowsamp\] is that, with high probability, a best $k$-rank approximation of a given matrix based on a subspace spanned by the sampled row is with in a relative $\epsilon$ error to the best $k$-rank approximation given by SVD. We list a few methods that use this observation. However, the complexity of finding a particular $k$-rank approximation to an $m\times n$ matrix is still ${\mathcal{O}}(kmn)$, as the complexity truncated SVD algorithms using Arnoldi or Lanczos methods [@GolV96; @LSY98]. In §\[Sub:CURap\] we recall the CUR-approximation introduced in [@GTZ97]. The main idea of CUR-approximation is to choose $k$ columns and rows of $A$, viewed as matrices $C$ and $R$, and then to choose a square matrix $U$ of order $k$ in such a way that $CUR$ is an optimal approximation of $A$. The matrix $U$ is chosen to be the inverse of the corresponding $k\times k$ submatrix $A'$ of $A$. The quality of CUR-approximation can be determined by the ratio of $|\det A'|$ to the maximum possible value of the absolute value of all $k\times k$ minors of $A$. In practice one searches for this maximum using a number of random choices of such minors. A modification of this search algorithm is given in [@FMMN11]. The complexity of storage of $C,R,U$ is ${\mathcal{O}}(k\max(m,n))$. The complexity of finding the value of each entry of $CUR$ is ${\mathcal{O}}(k^2)$. The complexity of computation of $CUR$ is ${\mathcal{O}}(k^2mn)$. In §\[S:fastapten\] we survey CUR-approximation of tensors given in [@FMMN11]. In §\[S:prelbrap\] we discuss preliminary results on best ${\mathbf{r}}$-approximation of tensors. In §\[sub:maxprob\] we show that the minimum problem is equivalent to the maximum problem . In §\[sub:singvten\] we discuss the notion of singular tuples and singular valuesva tensor introduced in [@Lim05]. In §\[sub:apprprb\] we recall the well known solution of maximizing $\|P_{\otimes_{i\in[d]}{\mathbf{U}}_i}({\mathcal{T}})\|^2$ with respect to one subspace, while keeping other subspaces fixed. In §\[S:AMM\] we discuss AMM for best ${\mathbf{r}}$-approximation and its variations. (In [@LMV00; @FMPS13] AMM is called *alternating least squares*, abbreviated as ALS.) In §\[sub:gendp\] we discuss the AMM on a product space. We mention a *modified alternating maximization method* and and *2-alternating maximization method*, abbreviated as MAMM and 2AMM respectively, introduced in [@FMPS13]. The MAMM method consists of choosing the one variable which gives the steepest ascend of AMM. 2AMM consists of maximization with respect to a pair of variables, while keeping all other variables fixed. In §\[sub:AMMbrapr\] we discuss briefly AMM and MAMM for best ${\mathbf{r}}$-approximations for tensors. In §\[sub:companbrap\] we give the complexity analysis of AMM for $d=3, r_1\approx r_2\approx r_3$ and $n_1\approx n_2\approx n_3$. In §\[S:fixpt\] we state a working assumption of this paper that AMM converges to a fixed point of the induced map, which satisfies certain smoothness assumptions. Under these assumptions we can apply the Newton method, which can be stated in the standard form in ${\mathbb{R}}^L$. Thus, we first do a number of AMM iterations and then switch to the Newton method. In §\[sub:newtboneappr\] we give a simple application of these ideas to state a Newton method for best rank one approximation. This Newton method was suggested in [@FMPS11]. It is different from the Newton method in [@ZhaG01]. The new contribution of this paper is the Newton method which is discussed in §\[S:bestbrapprox\] and §\[S:forDF\]. The advantage of our Newton method is its simplicity, which avoids the notions and tools of Riemannian geometry as for example in [@ES09; @SL10]. In simulations that we ran, the Newton method in [@ES09] was $20\%$ faster than our Newton method for best ${\mathbf{r}}$-approximation of $3$-mode tensors. However, the number of iterations of our Newton method was $40\%$ less than in [@ES09]. In the last section we give numerical results of our methods for best ${\mathbf{r}}$-approximation of tensors. In §\[S:NumRes\] we give numerical simulations of our different methods applied to a real computer tomography (CT) data set (the so-called MELANIX data set of OsiriX). The summary of these results are given in §\[S:Conclusion\].

Singular Value Decomposition {#S:SVD}
============================

Let $A\in{\mathbb{R}}^{m\times n}\setminus\{0\})$. We now recall well known facts on the SVD of $A$ [@GolV96]. See [@Ste93] for the early history of the SVD. Assume that $r={\mathrm{rank\;}}A$. Then there exist $r$-orthonormal sets of vectors ${\mathbf{u}}_1,\ldots,{\mathbf{u}}_r\in {\mathbb{R}}^m, {\mathbf{v}}_1,\ldots,{\mathbf{v}}_r\in {\mathbb{R}}^n$ such that we have: $$\begin{aligned}
\notag
&&A{\mathbf{v}}_i=\sigma_i(A){\mathbf{u}}_i, \quad {\mathbf{u}}_i{^\top}A=\sigma_i(A){\mathbf{v}}_i{^\top}, \quad i\in [r], \quad \sigma_1(A)\ge \cdots\ge \sigma_r(A)>0,\\
&&A_k=\sum_{i\in[k]} \sigma_i(A){\mathbf{u}}_i{\mathbf{v}}_i{^\top}, \quad k\in [r], \; A=A_r. \label{SVDid}\end{aligned}$$ The quantities ${\mathbf{u}}_i$, ${\mathbf{v}}_i$ and $\sigma_i(A)$ are called the left, right $i$-th singular vectors and $i$-th singular value of $A$ respectively, for $i\in[r]$. Note that ${\mathbf{u}}_k$ and ${\mathbf{v}}_k$ are uniquely defined up to $\pm 1$ if and only if $\sigma_{k-1}(A)>\sigma_k(A)>\sigma_{k+1}(A)$. Furthermore for $k\in [r-1]$ the matrix $A_k$ is uniquely defined if and only if $\sigma_k(A)>\sigma_{k+1}(A)$. Denote by ${\mathcal{R}}(m,n,k)\subset {\mathbb{R}}^{m\times n}$ the variety of all matrices of rank at most $k$. Then $A_k$ is a best rank-$k$ approximation of $A$: $$\min_{B\in{\mathcal{R}}(m,n,k)} \|A-B\|=\|A-A_k\|.$$ Let ${\mathbf{U}}\in {\mathop{\mathrm{Gr}}\nolimits}(p,{\mathbb{R}}^m), {\mathbf{V}}\in{\mathop{\mathrm{Gr}}\nolimits}(q,{\mathbb{R}}^n)$. We identify ${\mathbf{U}}\otimes {\mathbf{V}}$ with $$\label{defUVT}
{\mathbf{U}}{\mathbf{V}}{^\top}:={\mathrm{span}}\{{\mathbf{u}}{\mathbf{v}}{^\top},\quad {\mathbf{u}}\in{\mathbf{U}},{\mathbf{v}}\in {\mathbf{V}}\}\subset {\mathbb{R}}^{m\times n}.$$ Then $P_{{\mathbf{U}}\otimes {\mathbf{V}}}(A)$ is identified with the projection of $A$ on ${\mathbf{U}}{\mathbf{V}}{^\top}$ with respect to the standard inner product on ${\mathbb{R}}^{m\times n}$ given by ${\lanX,Y{\rangle}}={\mathop{\mathrm{tr}}\nolimits}XY{^\top}$. Observe that $$\textrm{Range }A={\mathbf{U}}_r^\star, \; {\mathbb{R}}^m={\mathbf{U}}_r^\star\oplus({\mathbf{U}}_r^\star)^\perp,  \textrm{ Range }A{^\top}={\mathbf{V}}_r^\star,\;
{\mathbb{R}}^n={\mathbf{V}}_r^\star\oplus({\mathbf{V}}_r^\star)^\perp.$$ Hence $$P_{{\mathbf{U}}\otimes {\mathbf{V}}}(A)=P_{({\mathbf{U}}\cap{\mathbf{U}}_r^\star)\otimes({\mathbf{V}}\cap{\mathbf{V}}_r^\star)}(A)\Rightarrow {\mathrm{rank\;}}P_{{\mathbf{U}}\otimes {\mathbf{V}}}(A)\le 
\min(\dim {\mathbf{U}},\dim {\mathbf{V}},r).$$ Thus $$\begin{aligned}
\notag
&&\max_{{\mathbf{U}}\in{\mathop{\mathrm{Gr}}\nolimits}(p,{\mathbb{R}}^m),{\mathbf{V}}\in{\mathop{\mathrm{Gr}}\nolimits}(q,{\mathbb{R}}^n)}\|P_{{\mathbf{U}}\otimes {\mathbf{V}}}(A)\|^2=\|P_{{\mathbf{U}}_l^\star\otimes{\mathbf{V}}_l^\star}(A)\|^2=\sum_{j\in[l]}\sigma_j(A)^2,\\
&&\min_{{\mathbf{U}}\in{\mathop{\mathrm{Gr}}\nolimits}(p,{\mathbb{R}}^m),{\mathbf{V}}\in{\mathop{\mathrm{Gr}}\nolimits}(q,{\mathbb{R}}^n)}\|A-P_{{\mathbf{U}}\otimes {\mathbf{V}}}(A)\|^2=\|A-P_{{\mathbf{U}}_l^\star\otimes{\mathbf{V}}_l^\star}(A)\|^2=\sum_{j\in[r]\setminus[l]}\sigma_j(A)^2,
\notag\\
&&l=\min(p,q,r). \label{maxprojmat}\end{aligned}$$

To compute ${\mathbf{U}}_l^\star,{\mathbf{V}}_l^\star$ and $\sigma_1(A),\ldots,\sigma_l(A)$ of a large scale matrix $A$ one can use Arnoldi or Lanczos methods [@GolV96; @LSY98], which are implemented in the partial singular value decomposition. This requires a substantial number of matrix-vector multiplications with the matrix $A$ and thus a complexity of at least ${\mathcal{O}}(lmn)$.

Sampling in low rank approximation of matrices {#S:lowrankapmat}
==============================================

Let $A=[a_{i,j}]_{i=j=1}^{m,n}\in{\mathbb{R}}^{m\times n}$ be given. Assume that ${\mathbf{b}}_1,\ldots,{\mathbf{b}}_m\in{\mathbb{R}}^n, {\mathbf{c}}_1,\ldots,{\mathbf{c}}_n\in{\mathbb{R}}^m$ are the columns of $A{^\top}$ and $A$ respectively. (${\mathbf{b}}_1{^\top},\ldots{\mathbf{b}}_m{^\top}$ are the rows of $A$.) Most of the known fast rank $k$-approximation are using sampling of rows or columns of $A$, or both.

Low rank approximations using sampling of rows {#Sub:rowsamp}
----------------------------------------------

Suppose that we sample a set $$\label{defsubI}
I:=\{1\le i_1<\ldots<i_s\le m\}\subset [m], \quad |I|=s,$$ of rows ${\mathbf{b}}_{i_1}{^\top},\ldots,{\mathbf{b}}_{i_s}{^\top}$, where $s\ge k$. Let ${\mathbf{W}}(I):={\mathrm{span}}({\mathbf{b}}_{i_1},\ldots,{\mathbf{b}}_{i_s})$. Then with high probability the projection of the first $i$-th right singular vectors ${\mathbf{v}}_i$ on ${\mathbf{W}}(I)$ is very close to ${\mathbf{v}}_i$ for $i\in[k]$, provided that $s\gg k$. In particular, [@DV06 Theorem 2] claims:

\[DesVemtheo\] (Deshpande-Vempala) Any $A\in{\mathbb{R}}^{m\times n}$ contains a subset $I$ of $s=\frac{4k}{\epsilon}+2k \log(k +1)$ rows such that there is a matrix $\tilde A_k$ of rank at most $k$ whose rows lie in ${\mathbf{W}}(I)$ and $$\|A-\tilde A_k\|^2\le (1+\epsilon)\|A-A_k\|^2.$$

To find a rank-$k$ approximation of $A$, one projects each row of $A$ on ${\mathbf{W}}(I)$ to obtain the matrix $P_{{\mathbb{R}}^m\otimes {\mathbf{W}}(I)}(A)$. Note that we can view $P_{{\mathbb{R}}^m\otimes {\mathbf{W}}(I)}(A)$ as an $m\times s'$ matrix, where $s'=\dim{\mathbf{W}}(I)\le s$. Then find a best rank $k$-approximation to $P_{{\mathbb{R}}^m\otimes {\mathbf{W}}(I)}(A)$, denoted as $P_{{\mathbb{R}}^m\otimes {\mathbf{W}}(I)}(A)_k$. Theorem \[DesVemtheo\] and the results of [@FKV04] yield that $$\|A-P_{{\mathbb{R}}^m\otimes {\mathbf{W}}(I)}(A)_k\|^2\le (1+\epsilon)\|A-A_k\|^2 + \eta \|A-P_{{\mathbb{R}}^m\otimes {\mathbf{W}}(I)}(A)\|^2.$$ Here $\eta$ is proportional to $\frac{k}{s}$, and can be decreased with more rounds of sampling. Note that the complexity of computing $P_{{\mathbb{R}}^m\otimes {\mathbf{W}}(I)}(A)_k$ is ${\mathcal{O}}(ks'm)$. The key weakness of this method is that to compute $P_{{\mathbb{R}}^m\otimes {\mathbf{W}}(I)}(A)$ one needs ${\mathcal{O}}(s'mn)$ operations. Indeed, after having computed an orthonormal basis of ${\mathbf{W}}(I)$, to compute the projection of each row of $A$ on ${\mathbf{W}}(I)$ one needs $s'n$ multiplications.

An approach for finding low rank approximations of $A$ using random sampling of rows or columns is given in Friedland-Kave-Niknejad-Zare [@FriKNZ06]. Start with a random choice of $I$ rows of $A$, where $|I|\ge k$ and $\dim {\mathbf{W}}(I)\ge k$. Find $P_{{\mathbb{R}}^m\otimes {\mathbf{W}}(I)}(A)$ and $B_1:=P_{{\mathbb{R}}^m\otimes {\mathbf{W}}(I)}(A)_k$ as above. Let ${\mathbf{U}}_{1}\in{\mathop{\mathrm{Gr}}\nolimits}(k,{\mathbb{R}}^m)$ be the subspace spanned by the first $k$ left singular vectors of $B_1$. Find $B_2=P_{{\mathbf{U}}_1\otimes {\mathbb{R}}^n}(A)$. Let ${\mathbf{V}}_2\in{\mathop{\mathrm{Gr}}\nolimits}(k,{\mathbb{R}}^n)$ correspond to the first $k$ right singular vectors of $B_2$. Continuing in this manner we obtain a sequence of rank $k$-approximations $B_1,B_2,\ldots$. It is shown in [@FriKNZ06] that $\|A-B_1\|\ge \|A-B_2\|\ge \ldots$ and $\|B_1\|\le \|B_2\|\le \ldots$. One stops the iterations when the relative improvement of the approximation falls below the specified threshold. Assume that $\sigma_k(A)>\sigma_{k+1}(A)$. Since best rank-$k$ approximation is a unique local minimum for the function $\|A-B\|, B\in {\mathcal{R}}(m,n,k)$ [@GolV96], it follows that in general the sequence $B_j, j\in{\mathbb{N}}$ converges to $A_k$. It is straightforward to show that this algorithm is the AMM for low rank approximations given in §\[sub:AMMbrapr\]. Again, the complexity of this method is ${\mathcal{O}}(kmn)$.

Other suggested methods as [@AM01; @DKM06; @RV07; @Sar06] seem to have the same complexity ${\mathcal{O}}(kmn)$, since they project each row of $A$ on some $k$-dimensional subspace of ${\mathbb{R}}^n$.

CUR-approximations {#Sub:CURap}
------------------

Let $$\label{defsubJ}
J:=\{1\le j_1<\ldots<j_t\le n\}\subset [n], \quad |J|=t,$$ and $I\subset [m]$ as in be given. Denote by $A[I,J]:=[a_{i_p,j_q}]_{p=q=1}^{s,t}\in {\mathbb{R}}^{s\times t}$. CUR-approximation is based on sampling simultaneously the set of $I$ rows and $J$ columns of $A$ and the approximation matrix to $A(I,J,U)$ given by $$\label{defCUR}
A(I,J,U):=CUR, \quad C:=A[[m],J], \; R:=A[I,[n]], \; U\in {\mathbb{R}}^{t\times s}.$$ Once the sets $I$ and $J$ are chosen the approximation $A(I,J,U)$ depends on the choice of $U$. Clearly the row and the column spaces of $A(I,J,U)$ are contained in the row and column spaces of $A[I,[n]]$ and $A[[m],J]$ respectively. Note that to store the approximation $A(I,J,U)$ we need to store the matrices $C$, $R$ and $U$. The number of these entries is $tm+sn+st$. So if $n,m$ are of order $10^5$ and $s,t$ are of order $10^2$ the storages of $C,R,U$ can be done in Random Access Memory (RAM), while the entries of $A$ are stored in external memory. To compute an entry of $A(I,J,U)$, which is an approximation of the corresponding entry of $A$, we need $st$ flops.

Let ${\mathbf{U}}$ and ${\mathbf{V}}$ be subspaces spanned by the columns of $A[[m],J]$ and $A[I,[n]]{^\top}$ respectively. Then $A(I,J,U)\in {\mathbf{U}}{\mathbf{V}}{^\top}$, see .

Clearly, a best CUR approximation is chosen by the least squares principle: $$\begin{aligned}
\notag
&&A(I,J,U^\star)=A([m],J)U^\star A(I,[n]),\\ 
&&U^\star=\arg\min\{\|A-A([m],J)UA(I,[n])\|,\; U\in{\mathbb{R}}^{|J|\times |I|}\}.\label{bestCURap}\end{aligned}$$ The results in [@FriT07] show that the least squares solution of is given by: $$\label{Ustarfor}
U^\star=A([m],J)^\dagger A A(I,[n])^\dagger.$$ Here $F^\dagger$ denotes the Moore-Penrose pseudoinverse of a matrix $F$. Note that $U^\star$ is unique if and only if $$\label{uniqUsol}
{\mathrm{rank\;}}A[[m],J]=|J|, \quad {\mathrm{rank\;}}A[I,[n]]=|I|.$$ The complexity of computation of $A([m],J)^\dagger$ and $A(I,[n])^\dagger$ are ${\mathcal{O}}(t^2m)$ and ${\mathcal{O}}(s^2n)$ respectively. Because of the multiplication formula for $U^\star$, the complexity of computation of $U^\star$ is ${\mathcal{O}}(stmn)$.

One can significantly improve the computation of $U$, if one tries to best fit the entires of the submatrix $A[I',J']$ for given subsets $I'\subset [m], J'\subset [n]$. That is, let $$\begin{aligned}
\notag
&&U^\star(I',J'):=\arg\min\{\|A[I',J']-A(I',J)UA(I,J')\|,\; U\in{\mathbb{R}}^{|J|\times |I|}\}=\\
&&A[I',J]^\dagger A[I',J'] A^\dagger [I,J'].
\label{UstarforI'J'}\end{aligned}$$ (The last equality follows from .) The complexity of computation of $U^\star(I',J')$ is ${\mathcal{O}}(st|I'||J'|)$.

Suppose finally, that $I'=I$ and $J'$. Then and the properties of the Moore-Penrose inverse yield that $$\label{UstarforIJ}
U^\star(I,J)=A[I,J]^\dagger, \quad A(I,J,U^\star(I,J))= B(I,J):=A[[m],J]A[I,J]^\dagger A[I,[n]].$$ In particular $B(I,J)[I,J]=A[I,J]$. Hence $$\begin{aligned}
\notag
&&A[[m],J]=B(I,J)[[m],J],\; A[I,[n]]=B(I,J)[I,[n]] \textrm{ if } |I|=|J|=k \textrm{ and } \det A[I,J]\ne 0,\\
&&B(I,J)=A[[m],J]A[I,J]^{-1} A[I,[n]].\label{invertAIJ}\end{aligned}$$ The original $CUR$ approximation of rank $k$ has the form $B(I,J)$ given by [@GTZ97].

Assume that ${\mathrm{rank\;}}A\ge k$. We want to choose an approximation $B(I,J)$ of the form which gives a good approximation to $A$. It is possible to give an upper estimate for the maximum of the absolute values of the entries of $A-B(I,J)$ in terms of $\sigma_{k+1}(A)$, provided that $\det A[I,J]$ is relatively close to $$\label{4.13.defmuk}
 \mu_k:=\max_{I\subset [m],J\subset[n],|I|=|J|=} |\det A[I,J]|>0.$$ Let $$\label{4.13.defentwnrm}
 \|F\|_{\infty,e}:=\max_{i\in[m],j\in[n]} |f_{i,j}|, \quad F=[f_{i,j}]\in{\mathbb{R}}^{m\times n}.$$

The results of [@GTZ97; @GorT01] yield: $$\label{4.13.basin}
 \|A- B(I,J)\|_{\infty,e}\le
 \frac{(k+1)\mu_k}{\det A[I,J]}\sigma_{p+1}(A).$$ (See also [@Fri15 Chapter 4, §13].)

To find $\mu_k$ is probably an NP-hard problem in general [@Fri13a]. A standard way to find $\mu_k$ is either a random search or greedy search [@GOSTZ10; @Fri13a]. In the special case when $A$ is a symmetric positive definite matrix one can give the exact conditions when the greedy search gives a relatively good result [@Fri13a].

In the paper by Friedland-Mehrmann-Miedlar-Nkengla [@FMMN11] a good approximation $B(I,J)$ of the form is obtained by a random search on the maximum value of the product of the significant singular values of $A[I,J]$. The approximations found in this way are experimentally better than the approximations found by searching for $\mu_k$.

Fast approximation of tensors {#S:fastapten}
=============================

The fast approximation of tensors can be based on several decompositions of tensors such as: Tucker decomposition [@Tuc66]; matricizations of tensors, as unfolding and applying SVD one time or several time recursively, (see below); higher order singular value decomposition (HOSVD) [@deLdV00], Tensor-Train decompositions [@Ose09; @Ose11]; hierarchical Tucker decomposition [@Gra10; @Hac11]. A very recent survey [@GKT13] gives an overview on this dynamic field. In this paper we will discuss only the CUR-approximation.

CUR-approximations of tensors {#sub:CURap}
-----------------------------

Let ${\mathcal{T}}\in{\mathbb{R}}^{n_1\times \ldots n_d}$. In this subsection we denote the entries of ${\mathcal{T}}$ as ${\mathcal{T}}(i_1,\ldots,i_d)$ for $i_j\in [n_j]$ and $j\in [d]$. CUR-approximation of tensors is based on matricizations of tensors. The unfolding of ${\mathcal{T}}$ in the mode $l\in [d]$ consists of rearranging the entries of ${\mathcal{T}}$ as a matrix $T_l({\mathcal{T}})\in {\mathbb{R}}^{n_l\times N_l}$, where $N_l=\frac{\prod_{i\in [d]}n_i}{n_l}$. More general, let $K\cup L=[d]$ be a partition of $[d]$ into two disjoint nonempty sets. Denote $N(K)=\prod_{i\in K} n_i, N(L)=\prod_{j\in L} n_j$. Then unfolding ${\mathcal{T}}$ into the two modes $K$ and $L$ consists of rearranging the entires of ${\mathcal{T}}$ as a matrix $T(K,L,{\mathcal{T}})\in {\mathbb{R}}^{N(K)\times N(L)}$.

We now describe briefly the $CUR$-approximation of $3$ and $4$-tensors as described by Friedland-Mehrmann-Miedlar-Nkengla [@FMMN11]. (See [@MMD06] for another approach to

CUR-approximations for tensors.) We start with the case $d=3$. Let $I_i$ be a nonempty subset of $[n_i]$ for $i\in [3]$. Assume that the following conditions hold: $$|I_1|=k^2, \quad |I_2|=|I_3|=k, \quad J:=I_2\times I_3\subset [n_2]\times [n_3].$$ We identify $ [n_2]\times [n_3]$ with $[n_2n_3]$ using a lexicographical order. We now take the CUR-approximation of $T_1({\mathcal{T}})$ as given in : $$B(I_1,J)=T_1({\mathcal{T}})[[n_1],J]T_1({\mathcal{T}})[I_1,J]^{-1}T_1({\mathcal{T}})[I_1,[n_2n_3]].$$ We view $T_1({\mathcal{T}})[[n_1],J]$ as an $n_1\times k^2$ matrix. For each $\alpha_1\in I_1$ we view $T_1({\mathcal{T}})[\{\alpha_1\}, [n_2n_3]]$ as an $n_2\times n_3$ matrix $Q(\alpha_1):=[{\mathcal{T}}(\alpha_1,i_2,i_3)]_{i_2\in[n_2], i_3\in[n_3]}$. Let $R(\alpha_1)$ be the $CUR$-approximation of $Q(\alpha_1)$ based on the sets $I_2,I_3$: $$R(\alpha_1):=Q(\alpha_1)[[n_2],I_3] Q(\alpha_1)[I_2,I_3]^{-1} Q(\alpha_1)[I_2,[n_3]].$$ Let $F:=T_1({\mathcal{T}})[I_1,J]^{-1}\in {\mathbb{R}}^{k^2\times k^2}$. We view the entries of this matrix indexed by the row $(\alpha_2,\alpha_3)\in I_2\times I_3$ and column $\alpha_1\in I_1$. We write these entries as ${\mathcal{F}}(\alpha_1,\alpha_2,\alpha_3), \alpha_j\in I_j, j\in[3]$, which represent a tensor ${\mathcal{F}}\in {\mathbb{R}}^{I_1\times I_2 \times I_3}$. The entries of $Q(\alpha_1)[I_2,I_3]^{-1}$ are indexed by the row $\alpha_3\in I_3$ and column $\alpha_2\in I_2$. We write these entries as ${\mathcal{G}}(\alpha_1,\alpha_2,\alpha_3), \alpha_2\in I_2,\alpha_3\in I_3$, which represent a tensor ${\mathcal{G}}\in {\mathbb{R}}^{I_1\times I_2 \times I_3}$. Then the approximation tensor ${\mathcal{B}}=[{\mathcal{B}}(j_1,j_2,j_3)]\in {\mathbb{R}}^{n_1\times n_2\times n_3}$ is given by: $$\begin{aligned}
&&{\mathcal{B}}(i_1,i_2,i_3)=\sum_{\alpha_1\in I_1,\alpha_j,\beta_j\in I_j, j=2,3}\\ 
&&{\mathcal{T}}(i_1,\alpha_2,\alpha_3){\mathcal{F}}(\alpha_1,\alpha_2,\alpha_3)
{\mathcal{T}}(\alpha_1,j_2,\beta_3){\mathcal{G}}(\alpha_1,\beta_2,\beta_3){\mathcal{T}}(\alpha_1,\beta_2,j_3).\end{aligned}$$

We now discuss a CUR-approximation for $4$-tensors, i.e. $d=4$. Let ${\mathcal{T}}\in {\mathbb{R}}^{n_1\times n_2\times n_3\times n_4}$ and $K=\{1,2\},L=\{3,4\}$. The rows and columns of $X:=T(K,L,{\mathcal{T}})\in {\mathbb{R}}^{(n_1n_2)\times (n_3n_4)}$ are indexed by pairs $(i_1,i_2)$ and $(i_3,i_4)$ respectively. Let $$I_j\subset [n_j], \;  |I_j|=k,\; j\in [4], \quad  J_1:=I_1\times I_2,\;J_2:=I_3\times I_4.$$ First consider the CUR-approximation $X[[n_1n_2],J_2] X[J_1,J_2]^{-1} X[J_1,[n_3n_4]]$ viewed as tensor ${\mathcal{C}}\in{\mathbb{R}}^{n_1\times n_2\times n_3\times n_4}$. Denote by ${\mathcal{H}}(\alpha_1,\alpha_2,\alpha_3,\alpha_4)$ the $((\alpha_3,\alpha_4), (\alpha_1,\alpha_2))$ entry of the matrix $X[J_1,J_2]^{-1}$. So ${\mathcal{H}}\in {\mathbb{R}}^{I_1\times I_2\times I_3\times I_4}$. Then $${\mathcal{C}}(i_1,i_2,i_3,i_4)=\sum_{\alpha_j\in I_j,j\in[4]} {\mathcal{T}}(i_1,i_2,\alpha_3,\alpha_4){\mathcal{H}}(\alpha_1,\alpha_2,\alpha_3,\alpha_4){\mathcal{T}}(\alpha_1,\alpha_2,i_3,i_4).$$

For $\alpha_j\in I_j,j\in[4]$ view vectors $X[[n_1n_2],(\alpha_3,\alpha_4)]$ and $X[(\alpha_1,\alpha_2),[n_3n_4]]$ as matrices $Y(\alpha_3,\alpha_4)\in {\mathbb{R}}^{n_1\times n_2}$ and $Z(\alpha_1,\alpha_2)\in {\mathbb{R}}^{n_3\times n_4}$ respectively. Next we find the CUR-approximations to these two matrices using the subsets $(I_1,I_2)$ and $(I_3,I_4)$ respectively: $$\begin{aligned}
&&Y(\alpha_3,\alpha_4)[[n_1],I_2]Y(\alpha_3,\alpha_4)[I_1,I_2]^{-1}Y(\alpha_3,\alpha_4)[I_1,[n_2]],\\
&&Z(\alpha_1,\alpha_2)[[n_3],I_4]Z(\alpha_1,\alpha_2)[I_3,I_4]^{-1}Z(\alpha_1,\alpha_2)[I_3,[n_4]].\end{aligned}$$ We denote the entries of $Y(\alpha_3,\alpha_4)[I_1,I_2]^{-1}$ and $Z(\alpha_1,\alpha_2)[I_3,I_4]^{-1}$ by ${\mathcal{F}}(\alpha_1,\alpha_2,\alpha_3,\alpha_4), \alpha_1\in I_1,\alpha_2\in I_2$ and ${\mathcal{G}}(\alpha_1,\alpha_2,\alpha_3,\alpha_4), \alpha_3\in I_3,\alpha_4\in I_4$ respectively. Then the CUR-approximation tensor ${\mathcal{B}}$ of ${\mathcal{T}}$ is given by: $$\begin{aligned}
&&{\mathcal{B}}(i_1,i_2,i_3,i_4)=\sum_{\alpha_j,\beta_j\in I_j, j\in [4]} {\mathcal{T}}(i_1,\beta_2,\alpha_3,\alpha_4){\mathcal{F}}(\beta_1,\beta_2,\alpha_3,\alpha_4)
{\mathcal{T}}(\beta_1,i_2,\alpha_3,\alpha_4)\\
&&{\mathcal{H}}(\alpha_1,\alpha_2,\alpha_3,\alpha_4){\mathcal{T}}(\alpha_1,\alpha_2,i_3,\beta_4){\mathcal{G}}(\alpha_1,\alpha_2,\beta_3,\beta_4){\mathcal{T}}(\alpha_1,\alpha_2,\beta_3,i_4).\end{aligned}$$

We now discuss briefly the complexity of the storage and computing an entry of the CUR-approximation ${\mathcal{B}}$. Assume first that $d=3$. Then we need to store $k^2$ columns of the matrices $T_1({\mathcal{T}})$, $k^3$ columns of $T_2({\mathcal{T}})$ and $T_3({\mathcal{T}})$, and $k^4$ entries of the tensors ${\mathcal{F}}$ and ${\mathcal{G}}$. The total storage space is $k^2 n_1+k^3(n_2+n_3)+2k^4$. To compute each entry of ${\mathcal{B}}$ we need to perform $4k^6$ multiplications and $k^6$ additions.

Assume now that $d=4$. Then we need to store $k^3$ columns of $T_l({\mathcal{T}}), l\in[4]$ and $k^4$ entries of ${\mathcal{F}},{\mathcal{G}},{\mathcal{H}}$. Total storage needed is $k^3(n_1+n_2+n_3+n_4+3k)$. To compute each entry of ${\mathcal{B}}$ we need to perform $6k^8$ multiplications and $k^8$ additions.

Preliminary results on best ${\mathbf{r}}$-approximation {#S:prelbrap}
========================================================

The maximization problem {#sub:maxprob}
------------------------

We first show that the best approximation problem is equivalent to the maximum problem , see [@LMV00] and [@Hac12 §10.3]. The Pythagoras theorem yields that $$\|{\mathcal{T}}\|^2 = \|P_{\otimes_{i=1}^d {\mathbf{U}}_i}({\mathcal{T}})\|^2 + \|P_{(\otimes_{i=1}^d {\mathbf{U}}_i)^\perp}({\mathcal{T}})\|^2,\quad
 \|{\mathcal{T}}-P_{\otimes_{i=1}^d {\mathbf{U}}_i}({\mathcal{T}})\|^2= \|P_{(\otimes_{i=1}^d {\mathbf{U}}_i)^\perp}({\mathcal{T}})\|^2.$$ (Here $(\otimes_{i=1}^d {\mathbf{U}}_i)^\perp$ is the orthogonal complement of $\otimes_{i=1}^d {\mathbf{U}}_i$ in $\otimes_{i=1}^d {\mathbb{R}}^{n_i}$.) Hence $$\label{bestapproxid}
 \min_{{\mathbf{U}}_i\in{\mathop{\mathrm{Gr}}\nolimits}(r_i,{\mathbb{R}}^{n_i}), i\in[d]}\|{\mathcal{T}}-P_{\otimes_{i=1}^d {\mathbf{U}}_i}({\mathcal{T}})\|^2=\|{\mathcal{T}}\|^2-
\max_{{\mathbf{U}}_i\in {\mathop{\mathrm{Gr}}\nolimits}(r_i,{\mathbb{R}}^{n_i}), i\in [d]} \|P_{\otimes_{i=1}{\mathbf{U}}_i}({\mathcal{T}})\|^2.$$ This shows the equivalence of and .

Singular values and singular tuples of tensors {#sub:singvten}
----------------------------------------------

Let ${\mathrm{S}}(n)=\{{\mathbf{x}}\in{\mathbb{R}}^n, \;\|{\mathbf{x}}\|=1\}$. Note that one dimensional subspace ${\mathbf{U}}\in{\mathop{\mathrm{Gr}}\nolimits}(1,{\mathbb{R}}^n)$ is ${\mathrm{span}}({\mathbf{u}})$, where ${\mathbf{u}}\in{\mathrm{S}}(n)$. Let ${\mathrm{S}}({\mathbf{n}}):={\mathrm{S}}(n_1)\times\cdots\times{\mathrm{S}}(n_d)$. Then best rank one approximation problem for ${\mathcal{T}}\in{\mathbb{R}}^{{\mathbf{n}}}$ is equivalent to finding $$\label{axprobbr1ap} 
 \|{\mathcal{T}}\|_{\sigma}:=\max_{({\mathbf{x}}_1,\ldots,{\mathbf{x}}_d)\in{\mathrm{S}}({\mathbf{n}})} {\mathcal{T}}\times (\otimes_{i\in[d]} {\mathbf{x}}_i).$$ Let $f_{\mathcal{T}}:{\mathbb{R}}^{{\mathbf{n}}}\to {\mathbb{R}}$ is given by $f_{{\mathcal{T}}}({\mathcal{X}})={{\langle}{\mathcal{X}},{\mathcal{T}}{\rangle}}$. Denote by ${\mathrm{S}}'({\mathbf{n}})\subset {\mathbb{R}}^{\mathbf{n}}$ all rank one tensors of the form $\otimes_{i\in[d]} {\mathbf{x}}_i$, where $({\mathbf{x}}_1,\ldots,{\mathbf{x}}_n)\in{\mathrm{S}}({\mathbf{n}})$. Let $f_{\mathcal{T}}({\mathbf{x}}_1,\ldots,{\mathbf{x}}_d):=f_{\mathcal{T}}(\otimes_{i\in[d]}{\mathbf{x}}_i)$. Then the critical points of $f_{\mathcal{T}}|{\mathrm{S}}'({\mathbf{n}})$ are given by the Lagrange multipliers formulas [@Lim05]: $$\label{defsingvalstup}
 {\mathcal{T}}\times (\otimes_{j\in[d]\setminus\{i\}} {\mathbf{u}}_j)=\lambda {\mathbf{u}}_i,\quad i\in [d],\quad ({\mathbf{u}}_1,\ldots,{\mathbf{u}}_d)\in{\mathrm{S}}({\mathbf{n}}).$$ One calls $\lambda$ and $({\mathbf{u}}_1,\ldots,{\mathbf{u}}_d)$ a singular value and singular tuple of ${\mathcal{T}}$. For $d=2$ these are the singular values and singular vectors of ${\mathcal{T}}$. The number of *complex* singular values of a generic ${\mathcal{T}}$ is given in [@FO14]. This number increases exponentially with $d$. For example for $n_1=\cdots=n_d=2$ the number of distinct singular values is $d!$. (The number of real singular values as given by is bounded by the numbers given in [@FO14].)

Consider first the maximization problem of $f_{\mathcal{T}}({\mathbf{x}}_1,\ldots,{\mathbf{x}}_d)$ over ${\mathrm{S}}({\mathbf{n}})$ where we vary ${\mathbf{x}}_i\in{\mathrm{S}}(n_i)$ and keep the other variables fixed. This problem is equivalent to the maximization of the linear form ${\mathbf{x}}_i{^\top}({\mathcal{T}}\times (\otimes_{j\in [d]\setminus\{i\}}{\mathbf{x}}_j))$. Note that if ${\mathcal{T}}\times (\otimes_{j\in [d]\setminus\{i\}}{\mathbf{x}}_j)\ne {\mathbf{0}}$ then this maximum is achieved for ${\mathbf{x}}_i=\frac{1}{\|{\mathcal{T}}\times (\otimes_{j\in [d]\setminus\{i\}}{\mathbf{x}}_j)\|}{\mathcal{T}}\times (\otimes_{j\in [d]\setminus\{i\}}{\mathbf{x}}_j)$.

Consider second the maximization problem of $f_{\mathcal{T}}({\mathbf{x}}_1,\ldots,{\mathbf{x}}_d)$ over ${\mathrm{S}}({\mathbf{n}})$ where we vary $({\mathbf{x}}_i, {\mathbf{x}}_j)\in{\mathrm{S}}(n_i)\times{\mathrm{S}}(n_j), 1\le i < j\le d$ and keep the other variables fixed. This problem is equivalent to finding the first singular value and the corresponding right and left singular vectors of the matrix ${\mathcal{T}}\times (\otimes_{k\in[d]\setminus\{i,j\}}{\mathbf{x}}_k)$. This can be done by using use Arnoldi or Lanczos methods [@GolV96; @LSY98]. The complexity of this method is ${\mathcal{O}}(n_in_j)$, given the matrix ${\mathcal{T}}\times (\otimes_{k\in[d]\setminus\{i,j\}}{\mathbf{x}}_k)$.

A basic maximization problem for best ${\mathbf{r}}$-approximation {#sub:apprprb}
------------------------------------------------------------------

Denote by ${\mathrm{S}}_n\subset {\mathbb{R}}^{n\times n}$ the space of real symmetric matrices. For $A\in {\mathrm{S}}_n$ denote by $\lambda_1(A)\ge \ldots\ge\lambda_n(A)$ the eigenvalues of $A$ arranged in a decreasing order and repeated according to their multiplicities. Let ${\mathrm{O}}(n,k)\subset {\mathbb{R}}^{n\times k}$ be the set of all $n\times k$ matrices $X$ with $k$ orthonormal columns, i.e. $X{^\top}X=I_k$, where $I_k$ is $k\times k$ identity matrix. We view $X\in {\mathbb{R}}^{n\times k}$ as composed of $k$-columns $[{\mathbf{x}}_1\ldots{\mathbf{x}}_k]$. The column space of $X\in {\mathrm{O}}(n,k)$ corresponds to a $k$-dimensional subspace ${\mathbf{U}}\subset{\mathbb{R}}^n$. Note that ${\mathbf{U}}\in{\mathop{\mathrm{Gr}}\nolimits}(k,{\mathbb{R}}^n)$ is spanned by the orthonormal columns of a matrix $Y\in O(n,k)$ if and only if $Y=XO$, for some $O\in {\mathrm{O}}(k,k)$.

For $A\in {\mathrm{S}}_n$ one has the Ky-Fan maximal characterization [@HJ88 Cor. 4.3.18] $$\label{KyFanchar}
 \max_{[{\mathbf{x}}_1\ldots{\mathbf{x}}_k]\in {\mathrm{O}}(n,k)}\sum_{i=1}^k {\mathbf{x}}_i{^\top}A{\mathbf{x}}_i=\sum_{i=1}^k \lambda_i(A).$$ Equality holds if and only if the column space of $X=[{\mathbf{x}}_1\ldots{\mathbf{x}}_k]$ is a subspace spanned by $k$ eigenvectors corresponding to $k$-largest eigenvalues of $A$.

We now reformulate the maximum problem in terms of orthonormal bases of ${\mathbf{U}}_i, i\in [d]$. Let ${\mathbf{u}}_{1,i},\ldots,{\mathbf{u}}_{n_i,i}$ be an orthonormal basis of ${\mathbf{U}}_i$ for $i\in [d]$. Then $\otimes_{i=1}^d {\mathbf{u}}_{j_i,i}, j_i\in [n_i],i\in[d]$ is an orthonormal basis of $\otimes_{i=1}^d {\mathbf{U}}_i$. Hence $$\|P_{\otimes_{i=1}^d{\mathbf{U}}_i}({\mathcal{T}})\|^2=\sum_{j_i\in [n_i],i\in[d]}{{\langle}{\mathcal{T}},\otimes_{i=1}^d {\mathbf{u}}_{j_i,i}{\rangle}}^2.$$ Hence is equivalent to $$\begin{aligned}
\label{maxnormprb1}
 &&\max_{[{\mathbf{u}}_{1,i}\ldots{\mathbf{u}}_{r_i,i}]\in{\mathrm{O}}(n_i,r_i),i\in[d]} \sum_{j_i\in [n_i], i\in[d]}{{\langle}{\mathcal{T}},\otimes_{i=1}^d {\mathbf{u}}_{j_i,i}{\rangle}}^2=\\
 &&\max_{{\mathbf{U}}_i\in{\mathop{\mathrm{Gr}}\nolimits}(r_i,{\mathbb{R}}^{n_i}),i\in[d]} \|P_{\otimes_{i=1}^d {\mathbf{U}}_i}({\mathcal{T}})\|^2. \notag
 \end{aligned}$$

A simpler problem is to find $$\begin{aligned}
\label{maxnormprbbas}
 &&\max_{[{\mathbf{u}}_{1,i}\ldots{\mathbf{u}}_{r_i,i}]\in{\mathrm{O}}(n_i,r_i)} \sum_{j_i\in [n_i], i\in[d]}{{\langle}{\mathcal{T}},\otimes_{i=1}^d {\mathbf{u}}_{j_i,i}{\rangle}}^2=\\
 &&\max_{{\mathbf{U}}_i\in{\mathop{\mathrm{Gr}}\nolimits}(r_i,{\mathbb{R}}^{n_i})} \|P_{\otimes_{i=1}^d {\mathbf{U}}_i}({\mathcal{T}})\|^2, \notag
 \end{aligned}$$ for a fixed $i\in [d]$. Let $$\begin{aligned}
\notag
 &&{\underline U}:=({\mathbf{U}}_1,\ldots,{\mathbf{U}}_d)\in{\mathop{\mathrm{Gr}}\nolimits}({\mathbf{r}},{\mathbf{n}}),\\
 &&{\mathop{\mathrm{Gr}}\nolimits}_i({\mathbf{r}},{\mathbf{n}}):={\mathop{\mathrm{Gr}}\nolimits}(r_1,n_1)\times\ldots\times
 {\mathop{\mathrm{Gr}}\nolimits}(r_{i-1},n_{i-1})\times{\mathop{\mathrm{Gr}}\nolimits}(r_{i+1},n_{i+1})\times\ldots\times{\mathop{\mathrm{Gr}}\nolimits}(r_d,n_d),\notag\\
 &&{\underline U}_i:=({\mathbf{U}}_1,\ldots,{\mathbf{U}}_{i-1},{\mathbf{U}}_{i+1},\ldots,{\mathbf{U}}_d)\in {\mathop{\mathrm{Gr}}\nolimits}_i({\mathbf{r}},{\mathbf{n}}),\notag\\
 &&A_i({\underline U}_i):=\sum_{j_l\in [r_l],l\in [d]\setminus\{i\}} 
 ({\mathcal{T}}\times \otimes_{k\in [d]\setminus \{i\}} {\mathbf{u}}_{j_k,k})({\mathcal{T}}\times \otimes_{k\in [d]\setminus \{i\}} {\mathbf{u}}_{j_k,k}){^\top}.
\label{defAUi}\end{aligned}$$

The maximization problem reduces to the maximum problem with $A=A_i({\underline U}_i)$. Note that each $A_i({\underline U}_i)$ is a positive semi-definite matrix. Hence $\sigma_j(A_i({\underline U}_i))=\lambda_j(A_i({\underline U}_i))$ for $j\in [n_i]$. Thus the complexity to find the first $r_i$ eigenvectors of $A_i({\underline U}_i)$ is ${\mathcal{O}}(r_in_i^2)$. Denote by ${\mathbf{U}}_i({\underline U}_i)\in{\mathop{\mathrm{Gr}}\nolimits}(r_i,{\mathbb{R}}^{n_i})$ a subspace spanned by the first $r_i$ eigenvectors of $A_i({\underline U}_i)$. Note that this subspace is unique if and only if $$\label{uniqcondUiAi}
\lambda_{r_i}(A_i({\underline U}_i))>\lambda_{r_i+1}(A_i({\underline U}_i)).$$ Finally, if ${\mathbf{r}}={\mathbf{1}}_d$ then each $A_i({\underline U}_i)$ is a rank one matrix. Hence ${\mathbf{U}}_i({\underline U}_i)={\mathrm{span}}({\mathcal{T}}\times \otimes_{k\in [d]\setminus \{i\}} {\mathbf{u}}_{1,k})$. For more details see [@FM08].

Alternating maximization methods for best ${\mathbf{r}}$-approximation {#S:AMM}
======================================================================

General definition and properties {#sub:gendp}
---------------------------------

Let $\Psi_i$ be a compact smooth manifold for $i\in[d]$. Define $$\Psi:=\Psi_1\times \cdots \times \Psi_d, \; 
\hat \Psi_i=(\Psi_1\times \cdots \times\Psi_{i-1}\times \Psi_{i+1}\times \cdots\times\Psi_d) \textrm{ for }i\in [d].$$ We denote by $\psi_i$, $\psi=(\psi_1,\ldots,\psi_d)$ and $\hat\psi_i=(\psi_1,\ldots,\psi_{i-1},\psi_{i+1},\ldots,\psi_d)$ the points in $\Psi_i$, $\Psi$ and $\hat \Psi_i$ respectively. Identify $\psi$ with $(\psi_i,\hat\psi_i)$ for each $i\in[d]$. Assume that $f:\Psi\to{\mathbb{R}}$ is a continuous function with continuous first and second partial derivatives. (In our applications it may happen that $f$ has discontinuities in first and second partial derivatives.) We want to find the maximum value of $f$ and a corresponding maximum point $\psi^\star$: $$\label{maxprobfPsi}
\max_{\psi\in\Psi} f(\psi)=f(\psi^\star).$$ Usually, this is a hard problem, where $f$ has many critical points and a number of these critical points are local maximum points. In some cases, as best ${\mathbf{r}}$ approximation to a given tensor ${\mathcal{T}}\in{\mathbb{R}}^{\mathbf{n}}$, we can solve the maximization problem with respect to one variable $\psi_i$ for any fixed $\hat\psi_i$: $$\label{maxprobPsii}
\max_{\psi_i\in \Psi_i} f((\psi_i,\hat\psi_i))=f((\psi_i^\star(\hat\psi_i),\hat\psi_i)),$$ for each $i\in[d]$.

Then the *alternating maximization method*, abbreviated as AMM, is as follows. Assume that we start with an initial point $\psi^{(0)}=(\psi^{(0)}_1,\ldots,\psi^{(0)}_d)=(\psi_1^{(0)},\hat \psi_1^{(0,1)})$. Then we consider the maximal problem for $i=1$ and $\hat\psi_1:=\hat \psi_1^{(0,1)}$. This maximum is achieved for $\psi_1^{(1)}:=\psi_1^\star(\hat\psi_1^{(0,1)})$. Assume that the coordinates $\psi_1^{(1)},\ldots,\psi_{j}^{(1)}$ are already defined for $j\in [d-1]$. Let $\hat \psi_{j+1}^{(0,j+1)}:=
(\psi_1^{(1)},\ldots,\psi_j^{(1)},\psi_{j+2}^{(0)},\ldots,\psi_d^{(0)})$. Then we consider the maximum problem for $i=j+1$ and $\hat\psi_{j+1}:=\hat \psi_{j+1}j^{(0,j+1)}$. This maximum is achieved for $\psi_{j+1}^{(1)}:=\psi_{j+1}^\star(\hat\psi_{j+1}^{(0,j+1)})$. Executing these $d$ iterations we obtain $\psi^{(1)}:=(\psi^{(1)}_1,\ldots,\psi^{(1)}_d)$. Note that we have a sequence of inequalities: $$f(\psi^{(0)})\le f(\psi_1^{(1)},\hat\psi_1^{(0,1)})\le f(\psi_2^{(1)},\hat\psi_2^{(0,2)})\le\cdots \le f(\psi_d^{(1)},\hat\psi_d^{(0,d)})=f(\psi^{(1)}).$$ Replace $\psi^{(0)}$ with $\psi^{(1)}$ and continue these iterations to obtain a sequence $\psi^{(l)}=(\psi_1^{(l)},\ldots,\psi_d^{(l)})$ for $l=0,\ldots,N$. Clearly, $$\label{nondecpropamm}
f(\psi^{(l-1)})\le f(\psi^{(l)})  \textrm{ for }l\in{\mathbb{N}}\Rightarrow \lim_{l\to\infty} f(\psi^{(l)})=M.$$ Usually, the sequence $\psi^{(l)}, l=0,\ldots,$ will converge to $1$-semi maximum point $\phi=(\phi_1,\ldots,\phi_d)\in\Psi$. That is, $f(\phi)=\max_{\psi_i\in\Psi} f((\psi_i,\hat\phi_i))$ for $i\in [d]$. Note that if $f$ is differentiable at $\phi$ then $\phi$ is a critical point of $f$. Assume that $f$ is twice differentiable at $\phi$. Then $\phi$ does not have to be a local maximum point [@FMPS13 Appendix].

The *modified alternating maximization method*, abbreviated as MAMM, is as follows. Assume that we start with an initial point $\psi^{(0)}=(\psi^{(0)}_1,\ldots,\psi^{(0)}_d)$. Let $\psi^{(0)}=(\psi_i^{(0)},\hat \psi_i^{(0)})$ for $i\in[d]$. Compute $f_{i,0}=
\max_{\psi_i\in\Psi_i} f((\psi_i,\hat\psi_i^{(0)}))$ for $i\in [d]$. Let $j_1\in\arg\max_{i\in[d]} f_{i,0}$. Then $\psi^{(1)}=(\psi_j^\star(\hat\psi_{j_1}^{(0)}),\hat\psi_{j_1}^{(0)})$ and $f_1=f_{1,{j_1}}=f(\psi^{(1)})$. Note that it takes $d$ iterations to compute $\psi^{(1)}$. Now replace $\psi^{(0)}$ with $\psi^{(1)}$ and compute $f_{i,1}=
\max_{\psi_i\in\Psi_i} f((\psi_i,\hat\psi_i^{(1)}))$ for $i\in[d]\setminus\{j_1\}$. Continue as above to find $\psi^{(l)}$ for $l=2,\ldots,N$. Note that for $l\ge 2$ it takes $d-1$ iterations to determine $\psi^{(l)}$. Clearly, holds. It is shown in [@FMPS13] that the limit $\phi$ of each convergent subsequence of the points $\psi^{(j)}$ is $1$-semi maximum point of $f$.

In certain very special cases, as for best rank one approximation, we can solve the maximization problem with respect to any pair of variables $\psi_i,\psi_j$ for $1\le i<j\le d$, where $d\ge 3$ and all other variables are fixed. Let $$\begin{aligned}
&&\hat \Psi_{i,j}:=\Psi_1\times \cdots \times\Psi_{i-1}\times \Psi_{i+1}\times\cdots\times \Psi_{j-1}\times \Psi_{j+1}\times\cdots\times \Psi_d,\\
&&\hat\psi_{i,j}=(\psi_1,\ldots,\psi_{i-1},\psi_{i+1},\ldots,\psi_{j-1},\psi_{j+1},\ldots,\psi_d)\in \hat \Psi_{i,j},\quad \psi_{i,j}=(\psi_i,\psi_j)\in \Psi_i\times \Psi_j.\end{aligned}$$ View $\psi=(\psi_1,\ldots,\psi_d)$ as $(\psi_{i,j},\hat\psi_{i,j})$ for each pair $1\le i<j\le d$. Then $$\label{maxprobPsiij}
\max_{\psi_{i,j}\in \Psi_i\times \Psi_j} f((\psi_{i,j},\hat\psi_{i,j}))=f((\psi_{i,j}^\star(\hat\psi_{i,j}),\hat\psi_{i,j})).$$ A point $\psi$ is called *2-semi maximum point* if the above maximum equals to $f(\psi)$ for each pair $1\le i<j\le d$.

The *$2$-alternating maximization method*, abbreviated here as 2AMM, is as follows. Assume that we start with an initial point $\psi^{(0)}=(\psi^{(0)}_1,\ldots,\psi^{(0)}_d)$. Suppose first that $d=3$. Then we consider the maximization problem for $i=2, j=3$ and $\hat\psi_{2,3}=\psi_1^{(0)}$. Let $(\psi_2^{(0,1)},\psi_3^{(0,1)})=\psi_{2,3}^\star(\psi_1^{(0)})$. Next let $i=1,j=3$ and $\psi_{1,3}=\psi_2^{(0,1)}$. Then $(\psi_1^{(0,2)},\psi_3^{(0,2)})=\psi_{1,3}^\star(\psi_2^{(0,1)})$. Next let $i=1,2$ and $\hat\psi_{1,2}=\psi_3^{(0,2)}$. Then $\psi^{(1)}=(\hat\psi_{1,2}^\star(\psi_3^{(0,2)}),\psi_3^{(0,2)})$. Continue these iterations to obtain $\psi^{(l)}$ for $l=2,\ldots$. Again, holds. Usually the sequence $\psi^{(l)},l\in{\mathbb{N}}$ will converge to a 2-semi maximum point $\phi$. For $d\ge 4$ the 2AMM can be defined appropriately see [@FMPS13].

A *modified $2$-alternating maximization method*, abbreviated here as M2AMM, is as follows. Start with an initial point $\psi^{(0)}=(\psi^{(0)}_1,\ldots,\psi^{(0)}_d)$ viewed as $(\psi_{i,j}^{(0)},\hat\psi_{i,j}^{(0)})$, for each pair $1\le i<j\le d$. Let $f_{i,j,0}:=\max_{\psi_{i,j}\in \Psi_i\times \Psi_j} f((\psi_{i,j},\hat\psi_{i,j}^{(0)}))$. Assume that

$(i_1,j_1)\in\arg\max_{1\le i <j\le d} f_{i,j,0}$. Then $\psi^{(1)}=(\psi_{i_1,j_1}^\star(\hat\psi_{i_1,j_1}^{(0)}),\hat\psi_{i_1,j_1}^{(0)})$. Let $f_{i_1,j_1,1}:=f(\psi^{(1)})$. Note that it takes $d \choose 2$ iterations to compute $\psi^{(1)}$. Now replace $\psi^{(0)}$ with $\psi^{(1)}$ and compute $f_{i,j,1}=
\max_{\psi_{i,j}\in\Psi_i\times \Psi_j} f((\psi_{i,j},\hat\psi_{i,j}^{(1)}))$ for all pairs $1\le i<j\le d$ except the pair $(i_1,j_1)$. Continue as above to find $\psi^{(l)}$ for $l=2,\ldots,N$. Note that for $l\ge 2$ it takes ${d \choose 2}-1$ iterations to determine $\psi^{(l)}$. Clearly, holds. It is shown in [@FMPS13] that the limit $\phi$ of each convergent subsequence of the points $\psi^{(j)}$ is $2$-semi maximum point of $f$.

AMM for best ${\mathbf{r}}$-aproximations of tensors {#sub:AMMbrapr}
----------------------------------------------------

Let ${\mathcal{T}}\in {\mathbb{R}}^{\mathbf{n}}$. For best rank one approximation one searches for the maximum of the function $f_{\mathcal{T}}={\mathcal{T}}\times (\otimes_{i\in [d]}{\mathbf{x}}_i)$ on ${\mathrm{S}}({\mathbf{n}})$, as in . For best ${\mathbf{r}}$-approximation one searches for the maximum of the function $f_{\mathcal{T}}=\|P_{\otimes_{i\in[d]}{\mathbf{U}}_i}\|^2$ on ${\mathop{\mathrm{Gr}}\nolimits}({\mathbf{r}},{\mathbf{n}})$, as in . A solution to the basic maximization problem with respect to one subspace ${\mathbf{U}}_i$ is given in §\[sub:apprprb\].

The AMM for best ${\mathbf{r}}$-approximation were studied first by de Lathauwer-Moor-Vandewalle [@LMV00]. The AMM is called in [@LMV00] *alternating least squares*, abbreviated as ALS. A crucial problem is the starting point of AMM. A high order SVD, abbreviated as HOSVD, for ${\mathcal{T}}$, see [@deLdV00], gives a good starting point for AMM. That is, let $T_l({\mathcal{T}})\in{\mathbb{R}}^{n_l\times N_l}$ be the unfolded matrix of ${\mathcal{T}}$ in the mode $l$, as in §\[sub:CURap\]. Then ${\mathbf{U}}_l$ is the subspace spanned by the first $l$-left singular vectors of $T_l({\mathcal{T}})$. The complexity of computing ${\mathbf{U}}_l$ is ${\mathcal{O}}(r_lN)$, where $N=\prod_{i\in[d]} n_i$. Hence for large $N$ the complexity of computing partial HOSVD is high. Another approach is to choose the starting subspaces at random, and repeat the AMM for several choices of random starting points.

MAMM for best rank one approximation of tensors was introduced by Friedland-Mehrmann-Pajarola-Suter in [@FMPS13] by the name *modified alternating least squares*, abbreviated as MALS. 2AMM for best rank one approximation was introduced in [@FMPS13] by the name *alternating SVD*, abbreviated as ASVD. It follows from the observation that $A:={\mathcal{T}}\times (\otimes_{l\in[d]\setminus\{i,j\}}{\mathbf{x}}_l)$ is an $n_i\times n_j$ matrix. Hence the maximum of the bilinear form ${\mathbf{x}}{^\top}A{\mathbf{y}}$ on ${\mathrm{S}}((n_i,n_j))$ is $\sigma_1(A)$. See §\[S:SVD\]. M2AMM was introduced in [@FMPS13] by the names MASVD.

We now introduce the following variant of 2AMM for best ${\mathbf{r}}$-rank approximation, called *$2$-alternating maximization method variant* and abbreviated as 2AMMV. Consider the maximization problem for a pair of variables as in . Since for ${\mathbf{r}}\ne {\mathbf{1}}_d$ we do not have a closed solution to this problem, we apply the AMM for two variables $\psi_i$ and $\psi_j$, while keeping $\hat \psi_{i,j}$ fixed. We then continue as in 2AMM method.

Complexity analysis of AMM for best ${\mathbf{r}}$-approximation {#sub:companbrap}
----------------------------------------------------------------

Let ${\underline U}=({\mathbf{U}}_1,\ldots,{\mathbf{U}}_d)$. Assume that $$\begin{aligned}
\notag
&&{\mathbf{U}}_{i}={\mathrm{span}}({\mathbf{u}}_{1,i},\ldots,{\mathbf{u}}_{r_i,n_i}), \quad {\mathbf{U}}_i^\perp={\mathrm{span}}({\mathbf{u}}_{r_i+1},\ldots, {\mathbf{u}}_{n_i,i}),\\
&&{\mathbf{u}}_{j,i}{^\top}{\mathbf{u}}_{k,i}=\delta_{j,k},\;j,k\in [n_i], \quad i\in [d].\label{basisUi}\end{aligned}$$ For each $i\in[d]$ compute the symmetric positive semi-definite matrix $A_i({\underline U}_i)$ given by . For simplicity of exposition we give the complexity analysis for $d=3$. To compute $A_1({\underline U}_1)$ we need first to compute the vectors ${\mathcal{T}}\times ({\mathbf{u}}_{j_2,2}\otimes {\mathbf{u}}_{j_3,3})$ for $j_2\in[r_2]$ and $j_3\in[r_3]$. Each computation of such a vector has complexity ${\mathcal{O}}(N)$, where $N=n_1n_2n_3$. The number of such vectors is $r_2r_3$. To form the matrix $A_i({\underline U}_i)$ we need ${\mathcal{O}}(r_2r_3n_1^2)$ flops. To find the first $r_1$ eigenvectors of $A_1({\underline U}_1)$ we need ${\mathcal{O}}(r_1n_1^2)$ flops. Assuming that $n_1,n_2,n_3\approx n$ and $r_1,r_2,r_3\approx r$ we deduce that we need $O(r^2 n^3)$ flops to find the first $r_1$ orthonormal eigenvectors of $A_1({\underline U}_1)$ which span ${\mathbf{U}}_1({\underline U}_1)$. Hence the complexity of finding orthonormal bases of ${\mathbf{U}}_1({\underline U}_1)$ is ${\mathcal{O}}(r^2n^3)$, which is the complexity of computing $A_1({\underline U}_1)$. Hence the complexity of each step of AMM for best ${\mathbf{r}}$-approximation, i.e. computing $\psi^{(l)}$, is ${\mathcal{O}}(r^2n^3)$.

It is possible to reduce the complexity of AMM for best ${\mathbf{r}}$-approximation is to ${\mathcal{O}}(r n^3)$ if we compute and store the matrices ${\mathcal{T}}\times {\mathbf{u}}_{j_1,1},{\mathcal{T}}\times {\mathbf{u}}_{j_2,2},{\mathcal{T}}\times {\mathbf{u}}_{j_3,3}$. See §\[sub:compNewtrap\].

We now analyze the complexity of AMM for rank one approximation. In this case we need only to compute the vector of the form ${\mathbf{v}}_i:={\mathcal{T}}\times(\otimes_{j\in[d]\setminus\{i\}}{\mathbf{u}}_j)$ for each $i\in[d]$, where ${\mathbf{U}}_j={\mathrm{span}}({\mathbf{u}}_j)$ for $j\in[i]$. The computation of each ${\mathbf{v}}_i$ needs ${\mathcal{O}}((d-2)N)$ flops, where $N=\prod_{j\in[d]}$. Hence each step of AMM for best rank one approximation is ${\mathcal{O}}(d(d-2)N)$. So for $d=3$ and $n_1\approx n_2\approx n_3$ the complexity is ${\mathcal{O}}(n^3)$, which is the same complexity as above with $r=1$.

Fixed points of AMM and Newton method {#S:fixpt}
=====================================

Consider the AMM as described in §\[sub:gendp\]. Assume that the sequence $\psi^{(l)}, l\in{\mathbb{N}}$ converges to a point $\phi\in\Psi$. Then $\phi$ is a fixed point of the map: $$\label{defFPsi}
\tilde{\mathbf{F}}:\Psi\to \Psi, \quad \tilde{\mathbf{F}}=(\tilde F_1,\ldots,\tilde F_d),\; \tilde F_i:\Psi\to \Psi_i,\; \tilde F_i(\psi)=\psi_i^\star(\hat\psi_i),\; \psi=(\psi_i,\hat\psi_i),\; i\in[d].$$ In general, the map $\tilde {\mathbf{F}}$ is a multivalued map, since the maximum given in may be achieved at a number of points denoted by $\psi_i^\star(\hat\psi_i)$. In what follows we assume:

\[fixpass\] The AMM converges to a fixed point $\phi$ of $\tilde {\mathbf{F}}$ i.e. $\tilde{\mathbf{F}}(\phi)=\phi$, such that the following conditions hold:

1.  There is a connected open neighborhood $O\subset \Psi$ such that $\tilde{\mathbf{F}}:O\to O$ is one valued map.

2.  $\tilde{\mathbf{F}}$ is a contraction on $O$ with respect to some norm on $O$.

3.  $\tilde{\mathbf{F}}\in {\mathrm{C}}^2(O)$, i.e. $\tilde{\mathbf{F}}$ has two continuous partial derivatives in $O$.

4.  $O$ is diffeomorphic to an open subset in ${\mathbb{R}}^L$. That is, there exists a smooth one-to-one map $H:O\to {\mathbb{R}}^L$ such that the Jacobian $D(H)$ is invertible at each point $\psi\in O$.

Assume that the conditions of Assumption \[fixpass\] hold. Then the map $\tilde {\mathbf{F}}: O\to O$ can be represented as $${\mathbf{F}}:O_1\to O_1, \quad {\mathbf{F}}=H\circ \tilde {\mathbf{F}}\circ H^{-1}, \quad O_1=H(O).$$ Hence to find a fixed point of $\tilde {\mathbf{F}}$ in $O$ it is enough to find a fixed point of ${\mathbf{F}}$ in $O_1$. A fixed point of ${\mathbf{F}}$ is a zero point of the system $$\label{zeroGx}
{\mathbf{G}}({\mathbf{x}})={\mathbf{0}}, \quad G({\mathbf{x}}):={\mathbf{x}}-{\mathbf{F}}({\mathbf{x}}).$$ To find a zero of ${\mathbf{G}}$ we use the standard Newton method.

In this paper we propose new Newton methods. We make a few iterations of AMM and switch to a Newton method assuming that the conditions of Assumption \[fixpass\] hold as explained above. A fixed point of the map $\tilde{\mathbf{F}}$ for best rank one approximation induces a fixed point of map ${\mathbf{F}}: {\mathbb{R}}^{\mathbf{n}}\to {\mathbb{R}}^{\mathbf{n}}$ [@FMPS13 Lemma 2]. Then the corresponding Newton method to find a zero of ${\mathbf{G}}$ is straightforward to state and implement, as explained in the next subsection. This Newton method was given in [@FMPS11 §5]. See also Zhang-Golub [@ZhaG01] for a different Newton method for best $(1,1,1)$ approximation.

Let $\tilde{\mathbf{F}}:{\mathop{\mathrm{Gr}}\nolimits}({\mathbf{r}},{\mathbf{n}})\to {\mathop{\mathrm{Gr}}\nolimits}({\mathbf{r}},{\mathbf{n}})$ be the induced map AMM. Each ${\mathop{\mathrm{Gr}}\nolimits}(r,{\mathbb{R}}^n)$ can be decomposed as a compact manifold to a finite number of charts as explained in §\[S:bestbrapprox\]. These charts induce standard charts of ${\mathop{\mathrm{Gr}}\nolimits}({\mathbf{r}},{\mathbf{n}})$. After a few AMM iterations we assume that the neighborhood $O$ of the fixed point of $\tilde{\mathbf{F}}$ lies in one the charts of ${\mathop{\mathrm{Gr}}\nolimits}({\mathbf{r}},{\mathbf{n}})$. We then construct the corresponding map ${\mathbf{F}}$ in this chart. Next we apply the standard Newton method to ${\mathbf{G}}$. The papers by Eldén-Savas [@ES09] and Savas-Lim [@SL10] discuss Newton and quasi-Newton methods for $(r_1,r_2,r_3)$ approximation of $3$-tensors using the concepts of differential geometry.

Newton method for best rank one approximation {#sub:newtboneappr}
---------------------------------------------

Let ${\mathcal{T}}\in{\mathbb{R}}^{{\mathbf{n}}}\setminus \{0\}$. Define: $$\begin{aligned}
\notag
&&\Psi_i={\mathbb{R}}^{n_i},\; i\in [d],\quad \Psi={\mathbb{R}}^{n_1}\times \cdots \times {\mathbb{R}}^{n_d}, \quad \psi=({\mathbf{x}}_1,\ldots,{\mathbf{x}}_d)\in \Psi,\\
&&f_{\mathcal{T}}:\Psi\to {\mathbb{R}}, \quad f_{\mathcal{T}}(\psi)={\mathcal{T}}\times (\otimes_{j\in[d]}{\mathbf{x}}_j),\\
&&{\mathbf{F}}=(F_1,\ldots,F_d):\Psi\to \Psi, \quad F_i(\psi)={\mathcal{T}}\times (\otimes_{j\in[d]\setminus\{i\}}{\mathbf{x}}_j), \quad i\in [d].\label{defbFfcT}\end{aligned}$$ Recall the results of §\[sub:singvten\]: Any critical point of $f_{\mathcal{T}}|{\mathrm{S}}({\mathbf{n}})$ satisfies . Suppose we start the AMM with $\psi^{(0)}=({\mathbf{x}}_1^{(0)},\ldots,{\mathbf{x}}^{(0)}_d)\in{\mathrm{S}}({\mathbf{n}})$ such that $f_{\mathcal{T}}(\psi^{(0)})\ne 0$. Then it is straightforward to see that $f_{\mathcal{T}}(\psi^{(l)})>0$ for $l\in{\mathbb{N}}$. Assume that $\lim_{l\to\infty}\psi^{(l)}=\omega=({\mathbf{u}}_1,\ldots,{\mathbf{u}}_d)
\in{\mathrm{S}}({\mathbf{n}})$. Then $\omega$ is the singular tuple of ${\mathcal{T}}$ satisfying . Clearly, $\lambda=f_{\mathcal{T}}(\omega)>0$. Let $$\label{defnorphi}
\phi=({\mathbf{y}}_1,\ldots,{\mathbf{y}}_d):=\lambda^{-\frac{1}{d-2}}\omega=\lambda^{-\frac{1}{d-2}}({\mathbf{u}}_1,\ldots,{\mathbf{u}}_d).$$ Then $\phi$ is a fixed point of ${\mathbf{F}}$.

Our Newton algorithm for finding the fixed point $\phi$ of ${\mathbf{F}}$ corresponding to a fixed point $\omega$ of AMM is as follows. We do a number of iterations of AMM to obtain $\psi^{(m)}$. Then we renormalize $\psi^{(m)}$ according to : $$\label{phi0for}
\phi_0:=(f_{\mathcal{T}}(\psi^{(m)})^{-\frac{1}{d-2}}\psi^{(m)}.$$ Let $D{\mathbf{F}}(\psi)$ denote the Jacobian of ${\mathbf{F}}$ at $\psi$, i.e. the matrix of partial derivatives of ${\mathbf{F}}$ at $\psi$. Then we perform Newton iterations of the form: $$\label{Newtitomeg}
\phi^{(l)}=\phi^{(l-1)} -(I-D{\mathbf{F}}(\phi^{(l-1)}))^{-1}(\phi^{(l-1)}-{\mathbf{F}}(\phi^{(l-1)})), \quad l\in{\mathbb{N}}.$$

After performing a number of Newton iterations we obtain $\phi^{(m')}=({\mathbf{z}}_1,\ldots,{\mathbf{z}}_d)$ which is an approximation of $\phi$. We then renormalize each ${\mathbf{z}}_i$ to obtain $\omega^{(m')}:=(\frac{1}{\|{\mathbf{z}}_1\|}{\mathbf{z}}_1,\ldots,\frac{1}{\|{\mathbf{z}}_d\|}{\mathbf{z}}_d)$ which is an approximation to the fixed point $\omega$. We call this Newton method *Newton-1*.

We now give the explicit formulas for $3$-tensors, where $n_1=m,n_2=n,n_3=l$. First $$\label{defmapF}
{\mathbf{F}}({\mathbf{u}},{\mathbf{v}},{\mathbf{w}}):=({\mathcal{T}}\times({\mathbf{v}}\otimes{\mathbf{w}}),{\mathcal{T}}\times({\mathbf{u}}\otimes{\mathbf{w}}),{\mathcal{T}}\times({\mathbf{u}}\otimes{\mathbf{v}})), \quad {\mathbf{G}}:=({\mathbf{u}},{\mathbf{v}},{\mathbf{w}})-{\mathbf{F}}({\mathbf{u}},{\mathbf{v}},{\mathbf{w}}).$$ Then $$\label{JacDF}
 D{\mathbf{G}}({\mathbf{u}},{\mathbf{v}},{\mathbf{w}})=\left[\begin{array}{ccc}I_m&-{\mathcal{T}}\times {\mathbf{w}}&-{\mathcal{T}}\times{\mathbf{v}}\\-({\mathcal{T}}\times{\mathbf{w}}){^\top}&I_n&-{\mathcal{T}}\times{\mathbf{u}}\\
 -({\mathcal{T}}\times {\mathbf{v}}){^\top}&-({\mathcal{T}}\times {\mathbf{u}}){^\top}&I_l\end{array}\right].$$ Hence Newton-1 iteration is given by the formula $$({\mathbf{u}}_{i+1},{\mathbf{v}}_{i+1},{\mathbf{w}}_{i+1})= ({\mathbf{u}}_{i},{\mathbf{v}}_{i},{\mathbf{w}}_i) -(D{\mathbf{G}}({\mathbf{u}}_i,{\mathbf{v}}_i,{\mathbf{w}}_i))^{-1}{\mathbf{G}}({\mathbf{u}}_i,{\mathbf{v}}_i,{\mathbf{w}}_i),$$ for $i=0,1,\ldots,$. Here we abuse notation by viewing $({\mathbf{u}},{\mathbf{v}},{\mathbf{w}})$ as a column vector $({\mathbf{u}}{^\top},{\mathbf{v}}{^\top},{\mathbf{w}}{^\top}){^\top}\in {\mathbb{C}}^{m+n+l}$.

Numerically, to find $(D{\mathbf{G}}({\mathbf{u}}_i,{\mathbf{v}}_i,{\mathbf{w}}_i))^{-1}{\mathbf{G}}({\mathbf{u}}_i,{\mathbf{v}}_i,{\mathbf{w}}_i)$ one solves the linear system $$(D{\mathbf{G}}({\mathbf{u}}_i,{\mathbf{v}}_i,{\mathbf{w}}_i))({\mathbf{x}},{\mathbf{y}},{\mathbf{z}})={\mathbf{G}}({\mathbf{u}}_i,{\mathbf{v}}_i,{\mathbf{w}}_i).$$ The final vector $({\mathbf{u}}_j,{\mathbf{v}}_j,{\mathbf{w}}_j)$ of Newton-1 iterations is followed by a scaling to vectors of unit length ${\mathbf{x}}_j=\frac{1}{\|{\mathbf{u}}_j\|}{\mathbf{u}}_j, {\mathbf{y}}_j=\frac{1}{\|{\mathbf{v}}_j\|}{\mathbf{v}}_j, {\mathbf{z}}_j=\frac{1}{\|{\mathbf{w}}_j\|}{\mathbf{w}}_j$.

We now discuss the complexity of Newton-1 method for $d=3$. Assuming that $m\approx n\approx l$ we deduce that the computation of the matrix $D{\mathbf{G}}$ is ${\mathcal{O}}(n^3)$. As the dimension of $D{\mathbf{G}}$ is $m+n+l$ it follows that the complexity of each iteration of Newton-1 method is ${\mathcal{O}}(n^3)$.

Newton method for best ${\mathbf{r}}$-approximation {#S:bestbrapprox}
===================================================

Recall that an $r$-dimensional subspace ${\mathbf{U}}\in {\mathop{\mathrm{Gr}}\nolimits}(r,{\mathbb{R}}^n)$ is given by a matrix $U=[u_{ij}]_{i,j=1}^{n,r}\in{\mathbb{R}}^{n\times r}$ of rank $r$. In particular there is a subset $\alpha\subset [n]$ of cardinality $r$ so that $\det U[\alpha,[r]]\ne 0$. Here $\alpha=(\alpha_1,\ldots,\alpha_r),
1\le \alpha_1<\ldots <\alpha_r\le n$. So $U[\alpha,[r]]:=[u_{\alpha_i j}]_{i,j=1}^r\in {\mathbf{GL}}(r,{\mathbb{R}})$, (the group of invertible matrices). Clearly, $V:=U U[\alpha,[r]]^{-1}$ represents another basis in ${\mathbf{U}}$. Note that $V[\alpha,[r]]=I_r$. Hence the set of all $V\in {\mathbb{R}}^{n\times r}$ with the condition: $V[\alpha,[r]]=I_r$ represent an open cell in ${\mathop{\mathrm{Gr}}\nolimits}(r,n)$ of dimension $r(n-r)$ denoted by ${\mathop{\mathrm{Gr}}\nolimits}(r,{\mathbb{R}}^n)(\alpha)$. (The number of free parameters in all such $V$’s is $(n-r)r$.) Assume for simplicity of exposition that $\alpha=[r]$. Note that $V_0=\left[\begin{array}{c}I_r\\0\end{array}\right]
\in {\mathop{\mathrm{Gr}}\nolimits}(r,{\mathbb{R}}^n)([r])$. Let ${\mathbf{e}}_i=(\delta_{1i},\ldots,\delta_{ni}){^\top}\in {\mathbb{R}}^n, i=1,\ldots,n$ be the standard basis in ${\mathbb{R}}^n$. So ${\mathbf{U}}_0={\mathrm{span}}({\mathbf{e}}_1,\ldots,{\mathbf{e}}_r)\in {\mathop{\mathrm{Gr}}\nolimits}(r,{\mathbb{R}}^n)([r])$, and $V_0$ is the unique representative of ${\mathbf{U}}_0$. Note that ${\mathbf{U}}_0^\perp$, the orthogonal complement of ${\mathbf{U}}_0$, is ${\mathrm{span}}({\mathbf{e}}_{r+1},\ldots,{\mathbf{e}}_n)$. It is straightforward to see that ${\mathbf{V}}\in {\mathop{\mathrm{Gr}}\nolimits}(r,{\mathbb{R}}^n)([r])$ if and only if ${\mathbf{V}}\cap {\mathbf{U}}_0^\perp=\{{\mathbf{0}}\}$.

The following definition is a geometric generalization of ${\mathop{\mathrm{Gr}}\nolimits}(r,{\mathbb{R}}^n)(\alpha)$: $$\label{openUcell}
{\mathop{\mathrm{Gr}}\nolimits}(r,{\mathbb{R}}^n)({\mathbf{U}}):=\{{\mathbf{V}}\in {\mathop{\mathrm{Gr}}\nolimits}(r,{\mathbb{R}}^n),\; {\mathbf{V}}\cap {\mathbf{U}}^\perp=\{{\mathbf{0}}\}\} \textrm{ for } {\mathbf{U}}\in{\mathop{\mathrm{Gr}}\nolimits}(r,{\mathbb{R}}^n).$$ A basis for ${\mathop{\mathrm{Gr}}\nolimits}(r,{\mathbb{R}}^n)({\mathbf{U}})$, which can be identified the tangent hyperplane $T_{\mathbf{U}}{\mathop{\mathrm{Gr}}\nolimits}(r,{\mathbb{R}}^n)$, can be represented as $\oplus^r {\mathbf{U}}^\perp$: Let ${\mathbf{u}}_1,\ldots,{\mathbf{u}}_r$ and ${\mathbf{u}}_{r+1},\ldots,{\mathbf{u}}_n$ be orthonormal bases of ${\mathbf{U}}$ and ${\mathbf{U}}^\perp$ respectively Then each subspace ${\mathbf{V}}\in {\mathop{\mathrm{Gr}}\nolimits}(r,{\mathbb{R}}^n)({\mathbf{U}})$ has a unique basis of the form ${\mathbf{u}}_1+{\mathbf{x}}_1,\ldots,{\mathbf{u}}_r+{\mathbf{x}}_r$ for unique ${\mathbf{x}}_1,\ldots,{\mathbf{x}}_r\in {\mathbf{U}}^\perp$. Equivalently, every matrix $X\in {\mathbb{R}}^{(n-r)\times r}$ induces a unique subspace ${\mathbf{V}}$ using the equality $$\label{defx1xr}
[{\mathbf{x}}_1\;\ldots\;{\mathbf{x}}_r]=[{\mathbf{u}}_1\;\ldots\;{\mathbf{u}}_{n-r}] X \textrm{ for each } X\in{\mathbb{R}}^{(n-r)\times r}.$$

Recall the results of §\[sub:apprprb\]. Let $\underline{U}=(U_1,\ldots,U_d)\in{\mathop{\mathrm{Gr}}\nolimits}({\mathbf{r}},{\mathbf{n}})$. Then $$\label{deftilFGrn}
\tilde{\mathbf{F}}=(\tilde F_1,\ldots, \tilde F_d):{\mathop{\mathrm{Gr}}\nolimits}({\mathbf{r}},{\mathbf{n}})\to{\mathop{\mathrm{Gr}}\nolimits}({\mathbf{r}},{\mathbf{n}}), \quad \tilde F_i({\underline U})={\mathbf{U}}_i({\underline U}_i), \; i\in[d],$$ where ${\mathbf{U}}_i({\underline U}_i)$ a subspace spanned by the first $r_i$ eigenvectors of $A_i({\underline U}_i)$. Assume that $\tilde{\mathbf{F}}$ is one valued at ${\underline U}$, i.e. holds. Then it is straightforward to show that $\tilde{\mathbf{F}}$ is smooth (real analytic) in neighborhood of $\underline{U}$. Assume next that there exists a neighborhood $O$ of $\underline{U}$ such that $$\label{regFcond}
 O\subset  {\mathop{\mathrm{Gr}}\nolimits}({\mathbf{r}},{\mathbf{n}})({\underline U}):={\mathop{\mathrm{Gr}}\nolimits}(r_1,{\mathbb{R}}^{n_1})({\mathbf{U}}_1)\times \cdots \times{\mathop{\mathrm{Gr}}\nolimits}(r_d,{\mathbb{R}}^{n_d})({\mathbf{U}}_d), \quad {\underline U}=({\mathbf{U}}_1,\ldots,{\mathbf{U}}_d),$$ such that the conditions *1-3* of Assumption \[fixpass\] hold. Observe next that ${\mathop{\mathrm{Gr}}\nolimits}({\mathbf{r}},{\mathbf{n}})({\underline U})$ is diffeomorphic to $${\mathbb{R}}^L:= {\mathbb{R}}^{(n_1-r_1)\times r_1} \ldots \times {\mathbb{R}}^{(n_d-r_d)\times r_d}, \quad L=\sum_{i\in[d]} (n_i-r_i)r_i$$ We say that $\tilde{\mathbf{F}}$ is *regular* at $\underline{U}$ if in addition to the above condition the matrix $I-D\tilde {\mathbf{F}}({\underline U})$ is invertible. We can view $X=[X_1\;\ldots \;X_d]\in  {\mathbb{R}}^{(n_1-r_1)\times r_1} \ldots \times {\mathbb{R}}^{(n_d-r_d)\times r_d}$. Then $\tilde {\mathbf{F}}$ on $O$ can be viewed as $$\label{tildebFre}
{\mathbf{F}}:O_1\to O_1, \quad O_1\subset {\mathbb{R}}^L,\quad {\mathbf{F}}(X)=[F_1(X),\ldots, F_d(X)],\; X=[X_1\;\ldots \;X_d]\in{\mathbb{R}}^L.$$ Note that $F_i(X)$ does not depend on $X_i$ for each $i\in[d]$. In our numerical simulations we first do a small number of AMM and then switch to Newton method given by . Observe that ${\underline U}$ corresponds to $X({\underline U})=[X_1({\underline U}),\ldots,X_d({\underline U})]$. When referring to we identify $X=[X_1,\ldots,X_d]$ with $\phi=(\phi_1,\ldots,\phi_d)$ and no ambiguity will arise.

Note that the case ${\mathbf{r}}={\mathbf{1}}_d$ corresponds to best rank one approximation. The above Newton method in this case is different from Newton method given in §\[sub:newtboneappr\].

A closed formula for $D{\mathbf{F}}(X({\underline U}))$ {#S:forDF}
=======================================================

Recall the definitions and results of §\[sub:apprprb\]. Given ${\underline U}$ we compute $\tilde F_i({\underline U})={\mathbf{U}}_i({\underline U}_i)$, which is the subspace spanned by the first $r_i$ eigenvectors of $A_i({\underline U}_i)$, which is given by , for $i\in[d]$. Assume that holds. Let $$\begin{aligned}
\notag
&&{\mathbf{U}}_{i}({\underline U}_i)={\mathrm{span}}({\mathbf{v}}_{1,i},\ldots,{\mathbf{v}}_{r_i,n_i}), \quad {\mathbf{U}}_i({\underline U}_i)^\perp={\mathrm{span}}({\mathbf{v}}_{r_i+1,i},\ldots, {\mathbf{v}}_{n_i,i}),\\
&&{\mathbf{v}}_{j,i}{^\top}{\mathbf{v}}_{k,i}=\delta_{jk},\;j,k\in [n_i], \quad i\in [d].\label{basistFUi}\end{aligned}$$

With each $X=[X_1,\ldots,X_d]\in{\mathbb{R}}^L$ we associate the following point $({\mathbf{W}}_1,\ldots,{\mathbf{W}}_d)\in{\mathop{\mathrm{Gr}}\nolimits}({\mathbf{r}},{\mathbf{n}})({\underline U})$. Suppose that $X_i=[x_{pq,i}]\in{\mathbb{R}}^{(n_i-r_i)\times r_i}$. Then ${\mathbf{W}}_i$ has a basis of the form $${\mathbf{u}}_{j_i,i}+\sum_{k_i\in [n_i-r_i]}x_{k_ij_i,i}{\mathbf{u}}_{r_i+k_i,i}, \quad j_i\in [r_i].$$ One can use the following notation for a basis ${\mathbf{w}}_{1,i},\ldots,{\mathbf{w}}_{r_i,i}$, written as a vector with vector coordinates $[{\mathbf{w}}_{1,i}\cdots{\mathbf{w}}_{r_i,i}]$: $$\label{basisWi}
[{\mathbf{w}}_{1,i}\cdots{\mathbf{w}}_{r_i,i}]=[{\mathbf{u}}_{1,i}\cdots{\mathbf{u}}_{r_i,i}] + [{\mathbf{u}}_{r_i+1,i}\cdots{\mathbf{u}}_{n_i,i}]X_i, \quad i\in[d].$$

Note that to the point ${\underline U}\in {\mathop{\mathrm{Gr}}\nolimits}({\mathbf{r}},{\mathbf{n}})({\underline U})$ corresponds the point $X=0$. Since ${\mathbf{u}}_{1,i},\ldots,{\mathbf{u}}_{n_i,i}$ is a basis in ${\mathbb{R}}^{n_i}$ it follows that $$\begin{aligned}
\notag
&&[{\mathbf{v}}_{1,i}\cdots {\mathbf{v}}_{r_i,i}]=[{\mathbf{u}}_{1,i}\cdots {\mathbf{u}}_{r_i,i}]Y_{i,0}+[{\mathbf{u}}_{r_i+1,i},\ldots,{\mathbf{u}}_{n_i,i}]X_{i,0}=[{\mathbf{u}}_{1,i}\cdots {\mathbf{u}}_{n_i,i}]Z_{i,0}, \\
&& Y_{i,0}\in {\mathbb{R}}^{r_i\times r_i}, \;X_{i,0}\in {\mathbb{R}}^{(n_i-r_i)\times r_i},\;Z_{i,0}=\left[\begin{array}{c}Y_{i,0}\\X_{i,0}\end{array}\right]\in {\mathbb{R}}^{n_i\times r_i},
\textrm{ for }i\in [d].\label{vforuXY}\end{aligned}$$ View $[{\mathbf{v}}_{1,i}\cdots {\mathbf{v}}_{r_i,i}],[{\mathbf{u}}_{1,i}\cdots {\mathbf{u}}_{n_i,i}]$ as $n_i\times r_i$ and $n_i\times n_i$ matrices with orthonormal columns. Then $$\label{forZi0}
Z_{i,0}=[{\mathbf{u}}_{1,i}\cdots {\mathbf{u}}_{n_i,i}]{^\top}[{\mathbf{v}}_{1,i}\cdots {\mathbf{v}}_{r_i,i}], \quad i\in [d].$$

The assumption that $\tilde {\mathbf{F}}:O\to O$ implies that $Y_{i,0}$ is an invertible matrix. Hence $[{\mathbf{v}}_{1,i}\cdots {\mathbf{v}}_{r_i,i}]Y_{i,0}^{-1}$ is also a basis in ${\mathbf{U}}_i({\underline U}_i)$. Clearly, $$[{\mathbf{v}}_{1,i}\cdots {\mathbf{v}}_{r_i,i}]Y_{i,0}^{-1}=[{\mathbf{u}}_{1,i}\cdots {\mathbf{u}}_{r_i,n_i}]+[{\mathbf{u}}_{r_i+1,i},\ldots,{\mathbf{u}}_{n_i,i}]X_{i,0}Y_{i,0}^{-1}, \quad i\in [d].$$

Hence $\tilde{\mathbf{F}}({\underline U})$ corresponds to ${\mathbf{F}}(0)$ where $$\label{exprebF0}
 F_i(0)=X_{i,0}Y_{i,0}^{-1}, \;i\in [d], \quad {\mathbf{F}}(0)=(F_1(0),\ldots,F_d(0)).$$

We now find the matrix of derivatives. So $D_iF_j\in{\mathbb{R}}^{((n_i-r_i)r_i\times ((n_j-r_j)r_j}$ is the partial derivative matrix of $(n_j-r_j)r_j$ coordinates of $F_j$ with respect to $(n_i-r_i)r_i$ the coordinates of ${\mathbf{U}}_i$ viewed as the matrix $\left[\begin{array}{c}I_{r_i}\\G_i\end{array}\right]$. So $G_i\in {\mathbb{R}}^{(n_i-r_i)\times r_i}$ are the variables representing the subspace ${\mathbf{U}}_i$. Observe first that $D_iF_i=0$ just as in Newton method for best rank one approximation in §\[sub:newtboneappr\].

Let us now find $D_i F_j(0)$. Recall that $D_i F_j(0)$ is a matrix of size $(n_i-r_i)r_i\times (n_j-r_j)r_j$. The entries of $D_i F_j(0)$ are indexed by $((p,q),(s,t))$ as follows: The entries of $G_i=[g_{pq,i}]\in {\mathbb{R}}^{(n_i-r_i)\times r_i}$ are viewed as $(n_i-r_i)r_i$ variables, and are indexed by $(p,q)$, where $p\in[n_i-r_i],q\in[r_i]$. $F_j$ is viewed as a matrix $G_j\in {\mathbb{R}}^{(n_j-r_j)\times r_j}$.The entries of $F_j$ are indexed by $(s,t)$, where $s\in[n_j-r_j]$ and $t\in[r_j]$. Since ${\underline U}\in{\mathop{\mathrm{Gr}}\nolimits}({\mathbf{r}},{\mathbf{n}})({\underline U})$ corresponds to $0\in{\mathbb{R}}^L$ we denote by $A_j(0)$ the matrix $A_j({\underline U}_j)$ for $j\in[d]$. We now give the formula for $\frac{\partial  A_j(0)}{\partial g_{pq,i}}$. This is done by noting that we vary ${\mathbf{U}}_{i}$ by changing the orthonormal basis ${\mathbf{u}}_{1,i},\ldots,{\mathbf{u}}_{r_i,i}$ up to the first perturbation with respect to the real variable $\varepsilon$ to $$\hat{\mathbf{u}}_{1,i}={\mathbf{u}}_{1,i},\ldots,\hat {\mathbf{u}}_{q-1,i}= {\mathbf{u}}_{q-1,i},\hat{\mathbf{u}}_{q,i}={\mathbf{u}}_{q,i}+\varepsilon{\mathbf{u}}_{r_i+p,i},
\hat{\mathbf{u}}_{q+1,i}={\mathbf{u}}_{q+1,i},\ldots,\hat{\mathbf{u}}_{r_i,i}={\mathbf{u}}_{r_i,i}$$ We denote the subspace spanned by these vectors as ${\mathbf{U}}_{i}(\varepsilon,p,q)$. That is, we change only the $q$ orthonormal vector of the standard basis in ${\mathbf{U}}_{i}$, for $q=1,\ldots,r_i$. The new basis is an orthogonal basis, and up order $\varepsilon$, the vector $ {\mathbf{u}}_{q,i}+\varepsilon{\mathbf{u}}_{r_i+p,i}$ is also of length $1$. Let ${\underline U}(\varepsilon,i,p,q)=({\mathbf{U}}_1,\ldots,{\mathbf{U}}_{i-1},{\mathbf{U}}_i(\varepsilon,p,q),{\mathbf{U}}_{i+1},\ldots,{\mathbf{U}}_d)$. Then ${\underline U}(\varepsilon,i,p,q)_j$ is obtained by dropping the subspace ${\mathbf{U}}_j$ from ${\underline U}(\varepsilon,i,p,q)$. We will show that $$\label{defAjvareps}
A_j({\underline U}(\varepsilon,i,p,q)_j)=A_j({\underline U}_j)+\varepsilon B_{j,i,p,q}+O(\varepsilon^2).$$

We now give a formula to compute $B_{j,i,p,q}$. Assume that $i,j\in [d]$ is a pair of different integers. Let $J$ be a set of $d-2$ pairs $\cup_{l\in [d]\setminus\{i,j\}}\{(k_l,l)\}$, where $k_l\in[r_l]$. Denote by ${\mathcal{J}}_{ij}$ the set of all such $J$’s. Note that ${\mathcal{J}}_{ij}={\mathcal{J}}_{ji}$. Furthermore, the number of elements in ${\mathcal{J}}_{ij}$ is $R_{ij}=\prod_{l\in[d]\setminus\{i,j\}} r_l$. We now introduce the following matrices $$\label{defCijk}
C_{ij}(J):={\mathcal{T}}\times (\otimes_{(k,l)\in J}{\mathbf{u}}_{k,l})\in {\mathbb{R}}^{n_i\times n_j}, \quad J\in{\mathcal{J}}_{ij}.$$ Note that $C_{ij}(J)=C_{ji}(J){^\top}$.

\[AjBjpqfor\] Let $i,j\in [d],i\ne j$. Assume that $p\in[n_i-r_i],q\in[r_i]$. Then holds. Furthermore $$\begin{aligned}
\label{newforAj}
&&A_j({\underline U}_j)=\sum_{k\in[r_i], J\in{\mathcal{J}}_{ji}} (C_{ji}(J){\mathbf{u}}_{k,i})(C_{ji}(J){\mathbf{u}}_{k,i}){^\top},\\
&&B_{j,i,p,q}=\sum_{J\in{\mathcal{J}}_{ji}} (C_{ji}(J){\mathbf{u}}_{k_i+p,i})(C_{ji}(J){\mathbf{u}}_{q,i}){^\top}+(C_{ji}(J){\mathbf{u}}_{q,i})(C_{ji}(J){\mathbf{u}}_{k_i+p,i}){^\top}\label{forBjipq}\end{aligned}$$

[**Proof. **]{}The identity of is just a restatement of . To compute $A_j({\underline U}(\varepsilon,i,p,q)_j)$ use by replacing $u_{k,i}$ with $\hat u_{k,i}$ for $k\in[n_i]$. Deduce first and then .[ $\Box$\
]{}

Recall that ${\mathbf{v}}_{1,j},\ldots,{\mathbf{v}}_{r_j,j}$ is an orthonormal basis of ${\mathbf{U}}_j({\underline U}_j)$, and these vectors are the eigenvectors $A_j({\underline U}_j)$ corresponding its first $r_j$ eigenvalues. Let ${\mathbf{v}}_{r_j+1,j},\ldots,{\mathbf{v}}_{n_j,i}$ be the last $n_j-r_j$ orthonormal eigenvectors of $A_j({\underline U}_j)$. We now find the first perturbation of the first $r_i$ eigenvectors for the matrix $A_j({\underline U}_j)+\varepsilon B_{j,i,p,q}$. Assume first, for simplicity of exposition, that each $\lambda_k(A_j({\underline U}_j))$ is simple for $k\in [r_j]$: Then it is known, e.g. [@Fri15 Chapter 4, §19, (4.19.2)]: $$\label{pertforxj2}
{\mathbf{v}}_{k,j}(\varepsilon,i,p,q)={\mathbf{v}}_{k,j}+\varepsilon (\lambda_{k}(A_j({\underline U}_j)) I_{n_j}-A_j({\underline U}_j))^\dagger B_{j,i,p,q}{\mathbf{v}}_{k,j} +O(\varepsilon^2),  k\in[r_j].$$

The assumption that $\lambda_{k}(A_j({\underline U}_j))$ is a simple eigenvalue for $k\in[r_j]$ yields $$(\lambda_{k}(A_j({\underline U}_i))I_{n_j}-A_j({\underline U}_j))^\dagger{\mathbf{y}}=\sum_{l\in[n_j] 
\setminus\{k\}}\frac{1}{\lambda_{k}(A_j({\underline U}_j))-\lambda_{l}(A_j({\underline U}_j))} ({\mathbf{v}}_{l,j}{^\top}{\mathbf{y}}){\mathbf{v}}_{l,j},$$ for ${\mathbf{y}}\in{\mathbb{R}}^{n_j}$.

Since we are interested in a basis of ${\mathbf{U}}_j({\underline U}(\varepsilon,i,p,q)_j)$ up to the order of $\varepsilon$ we can assume that this basis is of the form $$\tilde {\mathbf{v}}_{k,j}(\varepsilon,i,p,q)={\mathbf{v}}_{k,j}+\varepsilon{\mathbf{w}}_{k,j}(i,p,q), \quad {\mathbf{w}}_{k,j}(i,p,q)\in
{\mathrm{span}}({\mathbf{v}}_{r_j+1,j}\ldots,{\mathbf{v}}_{n_j,j}).$$ Hence $$\begin{aligned}
\notag
 &&{\mathbf{w}}_{k,j}(i,p,q)=\sum_{l\in[n_j]\setminus [r_j]}\frac{1}{\lambda_{k}(A_j({\underline U}_j))-\lambda_{l}(A_j({\underline U}_j))} ({\mathbf{v}}_{l,j}{^\top}{\mathbf{c}}_{k,j,i,p,q}){\mathbf{v}}_{l,j},\\
&&{\mathbf{c}}_{k,j,i,p,q}:= B_{j,i,p,q}{\mathbf{v}}_{k,j}.\label{pertforbasU2}\end{aligned}$$ Note that the assumption yields that ${\mathbf{w}}_{k,j}$ is well defined for $k\in[r_j]$. Let $$\begin{aligned}
&&W_j(i,p,q)=[{\mathbf{w}}_{1,j}(i,p,q)\cdots{\mathbf{w}}_{r_j,j}(i,p,q)]=\left[\begin{array}{c}V_j(i,p,q)\\U_j(i,p,q)\end{array}\right],\\ 
&&V_j(i,p,q)\in{\mathbb{R}}^{r_j\times r_j},\quad U_j(i,p,q)\in{\mathbb{R}}^{(n_j-r_j)\times r_j}.\end{aligned}$$ Up to the order of $\varepsilon$ we have that a basis of ${\mathbf{U}}_{j}({\underline U}(\varepsilon,i,p,q)_j)$ is given by columns of matrix $Z_{j,0}+\varepsilon W_j(i,p,q)=\left[\begin{array}{c}Y_{j,0}+\varepsilon V_j(i,p,q)\\X_{j,0} +\varepsilon U_j(i,p,q)\end{array}\right]$. Note $$(Z_{j,0}+\varepsilon W_j(i,p,q))(Y_{j,0}+\varepsilon V_j(i,p,q))^{-1}=\left[\begin{array}{c} I_{r_j}\\
(X_{j,0}+\varepsilon U_j(i,p,q))(Y_{j,0}+\varepsilon V_j(i,p,q))^{-1}\end{array}\right].$$ Observe next $$\begin{aligned}
&&Y_{j,0}+\varepsilon V_j(i,p,q)=Y_{j,0}(I_{r_j}+\varepsilon Y_{j,0}^{-1} V_j(i,p,q)), \\
&&(Y_{j,0}+\varepsilon V_j(i,p,q))^{-1}=(I_{r_j}+\varepsilon Y_{j,0}^{-1} V_j(i,p,q))^{-1}Y_{j,0}^{-1}=\\
&&Y_{j,0}^{-1}-\varepsilon Y_{j,0}^{-1} V_j(i,p,q)Y_{j,0}^{-1}+O(\varepsilon^2),\\
&&(X_{j,0}+\varepsilon U_j(i,p,q))(Y_{j,0}+\varepsilon V_j(i,p,q))^{-1}=\\
&&X_{j,0}Y_{j,0}^{-1}+\varepsilon(U_j(i,p,q) Y_{j,0}^{-1}-X_{j,0}Y_{j,0}^{-1}V_{j}(i,p,q) Y_{j,0}^{-1})+O(\varepsilon^2).\end{aligned}$$ Hence $$\label{D1F2pqfor}
\frac{\partial F_j}{\partial g_{pq,i}}(0)=U_j(i,p,q)Y_{j,0}^{-1}-X_{j,0}Y_{j,0}^{-1}V_{j}(i,p,q) Y_{j,0}^{-1}.$$

Thus $D{\mathbf{F}}(0)=[D_iF_J]_{i,j\in[d]}\in{\mathbb{R}}^{L\times L}$. We now make one iteration of Newton method given by for $l=1$, where $\phi^{(0)}=0$: $$\label{Newtmetstrt0}
\phi^{(1)}=-(I-D{\mathbf{F}}(0))^{-1}F(0), \quad \phi^{(1)}=[X_{1,1},\ldots,X_{d,1}]\in{\mathbb{R}}^L.$$ Let ${\mathbf{U}}_{i,1}\in{\mathop{\mathrm{Gr}}\nolimits}(r_i,{\mathbb{R}}^{n_i})$ be the subspace represented by the matrix $X_{i,1}$: $$\label{defUi1}
{\mathbf{U}}_{i,1}={\mathrm{span}}(\tilde{\mathbf{u}}_{1,i,1},\ldots,\tilde{\mathbf{u}}_{r_i,i,1}),\; [\tilde{\mathbf{u}}_{1,i,1},\ldots,\tilde{\mathbf{u}}_{n_i,i,1}]=
[{\mathbf{u}}_{1,i},\ldots,{\mathbf{u}}_{n_i,i}]\left[\begin{array}{c} I_{r_i}\\X_{j,1}\end{array}\right]$$ for $i\in[d]$. Perform the Gram-Schmidt process on $\tilde{\mathbf{u}}_{1,i,1},\ldots,\tilde{\mathbf{u}}_{r_i,i,1}$ to obtain an orthonormal basis ${\mathbf{u}}_{1,i,1},\ldots,{\mathbf{u}}_{r_i,i,1}$ of ${\mathbf{U}}_{i,1}$. Let ${\underline U}:=({\mathbf{U}}_{1,1},\ldots,{\mathbf{U}}_{d,1})$ and repeat the algorithm which is described above. We call this Newton method *Newton-2*.

Complexity of Newton-2 {#sub:compNewtrap}
======================

In this section we assume for simplicity that $d=3$, $r_1=r_2=r_3=r$, $n_i\approx n$ for $i\in [3]$. We assume that executed a number of times the AMM for a given ${\mathcal{T}}\in {\mathbb{R}}^{\mathbf{n}}$. So we are given ${\underline U}=({\mathbf{U}}_1,\ldots,{\mathbf{U}}_d)$, and an orthonormal basis ${\mathbf{u}}_{1,i},\ldots,{\mathbf{u}}_{r,i}$ of ${\mathbf{U}}_i$ for $i\in[d]$. First we complete each ${\mathbf{u}}_{1,i},\ldots,{\mathbf{u}}_{r,i}$ to an orthonormal basis ${\mathbf{u}}_{1,i},\ldots,{\mathbf{u}}_{n_i,i}$of ${\mathbb{R}}^{n_i}$, which needs ${\mathcal{O}}(n^3)$ flops. Since $d=3$ we still need only ${\mathcal{O}}(n^3)$ to carry out this completion for each $i\in [3]$.

Next we compute the matrices $C_{ij}(J)$. Since $d=3$, we need $n$ flops to compute each entry of $C_{ij}(J)$. Since we have roughly $n^2$ entries, the complexity of computing $C_{ij}(J)$ is ${\mathcal{O}}(n^3)$. As the cardinality of ${\mathcal{J}}_{ij}$ is $r$ we need ${\mathcal{O}}(rn^3)$ flops to compute all $C_{ij}(J)$ for $J\in{\mathcal{J}}_{ij}$. As the number of pairs in $[3]$ is $3$ it follows that the complexity of computing all $C_{ij}(J)$ is ${\mathcal{O}}(rn^3)$.

The identity yields that the complexity of computing $A_j({\underline U}_j)$ is ${\mathcal{O}}(r^2n^2)$. Recall next that $A_j({\underline U}_j)$ is $n_j\times n_j$ symmetric positive semi-definite matrix. The complexity of computations of the eigenvalues and the orthonormal eigenvectors of $A_j({\underline U}_j)$ is ${\mathcal{O}}(n^3)$. Hence the complexity of computing ${\underline U}$ is ${\mathcal{O}}(r n^3)$, as we pointed out at the end of §\[sub:companbrap\].

The complexity of computing $B_{j,i,p.q}$ using is ${\mathcal{O}}(r n^2)$. The complexity of computing ${\mathbf{w}}_{k,j}(i,p,q)$, given by is ${\mathcal{O}}(n^2)$. Hence the complexity of computing $W_j(i,p,q)$ is ${\mathcal{O}}(rn^2)$. Therefore the complexity of computing $D_i F_j$ is ${\mathcal{O}}(r^2n^3)$. Since $d=3$, the complexity of computing the matrix $D{\mathbf{F}}(0)$ is also ${\mathcal{O}}(r^2n^3)$.

As $D{\mathbf{F}}(0)\in {\mathbb{R}}^{L\times L}$, where $L\approx 3r n$, the complexity of computing $(I-D{\mathbf{F}}(0))^{-1}$ is ${\mathcal{O}}(r^3n^3)$. In summary, the complexity of one step in Newton-2 is ${\mathcal{O}}(r^3 n^3)$.

Numerical Results {#S:NumRes}
=================

We have implemented a Matlab library tensor decomposition using Tensor Toolbox given by [@KB09]. The performance was measured via the actual CPU-time (seconds) needed to compute. All performance tests have been carried out on a 2.8 GHz Quad-Core Intel Xeon Macintosh computer with 16GB RAM. The performance results are discussed for real data sets of third-order tensors. We worked with a real computer tomography (CT) data set (the so-called MELANIX data set of OsiriX) [@FMPS13].

Our simulation results are averaged over 10 different runs of the each algorithm. In each run, we changed the initial guess, that is, we generated new random start vectors. We always initialized the algorithms by random start vectors, because this is cheaper than the initialization via HOSVD. We note here that for Newton methods our initial guess is the subspaces returned by one iteration of AMM method.

All the alternating algorithms have the same stopping criterion where convergence is achieved if one of the two following conditions are met: $iterations > 10;fitchange <0.0001$ is met. All the Newton algorithms have the same stopping criterion where convergence is achieved if one of the two following conditions are met: $iterations > 10;change <\exp(-10)$.

Our numerical simulations demonstrate the well known fact that for large size tensors Newton methods are not efficient. Though the Newton methods converge in fewer iterations than alternating methods, the computation associated with the matrix of derivatives (Jacobian) in each iteration is too expensive making alternating maximization methods much more cost effective. Our simulations also demonstrate that our Newton-1 for best rank one approximation is as fast as AMM methods. However our Newton-2 is much slower than alternating methods. We also give a comparison between our Newton-2 and the Newton method based on Grassmannian manifold by [@ES09], abbreviated as Newton-ES.

We also observe that for large tensors and large rank approximation two alternating maximization methods, namely MAMM and 2AMMV, seem to outperform the other alternating maximization methods. We would recommend Newton-1 for rank one approximation in case of rank one approximation both for large and small sized tensors. For higher rank approximation we recommend 2AMMV in case of large size tensors and AMM or MAMM in case of small size tensors.

Our Newton-2 performs a bit slower than Newton-ES, however we would like to point couple of advantages. Our method can be easily extendable to higher dimensions ( for $d>3$ case) both analytically and numerically compared to Newton-ES. Our method is also highly parallelizable which can bring down the computation time drastically. Computation of $D_{i}F_{j}$ matrices in each iteration contributes to about $50\%$ of the total time, which however can be parallelizable. Finally the number of iterations in Newton-2 is at least $30\%$ less than in Newton-ES.

It is not only important to check how fast the different algorithms perform but also what quality they achieve. This was measured by checking the Hilbert-Schmidt norm, abbreviated as HS norm, of the resulting decompositions, which serves as a measure for the quality of the approximation. In general, we can say that the higher the HS norm, the more likely it is that we find a global maximum. Accordingly, we compared the HS norms to say whether the different algorithms converged to the same stationary point. In Figure 4, we show the average HS norms achieved by different algorithms and compared them with the input norm. We observe all the algorithms seem to attain the same local maximum.

Best (2,2,2) and rank two approximations {#sub:b222apr}
----------------------------------------

Assume that ${\mathcal{T}}$ is a $3$-tensor of rank three at least and let ${\mathcal{S}}$ be a best $(2,2,2)$-approximation to ${\mathcal{T}}$given by . It is easy to show that ${\mathcal{S}}$ has at least rank $2$. Let ${\mathcal{S}}'=[s_{j_1,j_2,j_3}]\in{\mathbb{R}}^{2\times 2 \times 2}$ be the core tensor corresponding to ${\mathcal{S}}$. Clearly ${\mathrm{rank\;}}{\mathcal{S}}={\mathrm{rank\;}}{\mathcal{S}}'\ge 2$. Recall that a real nonzero $2\times 2 \times 2$ tensor has rank one, two or three [@BK99]. So ${\mathrm{rank\;}}{\mathcal{S}}\in\{2,3\}$. Observe next that if ${\mathrm{rank\;}}{\mathcal{S}}={\mathrm{rank\;}}{\mathcal{S}}'=2$ then ${\mathcal{S}}$ is also a best rank two approximation of ${\mathcal{T}}$. Recall that a best rank two approximation of ${\mathcal{T}}$ may not always exist. In particular where ${\mathrm{rank\;}}{\mathcal{T}}>2$ and the border rank of ${\mathcal{T}}$ is $2$ [@DeSL08]. In all our numerical simulations for best $(2,2,2)$-approximation we performed on random large tensors, the tensor ${\mathcal{S}}'$ had rank two. Note that the probability of $2\times 2\times 2$ tensors, with entries normally distributed with mean $0$ and variance $1$, to have rank $2$ is $\frac{\pi}{4}$ [@Ber13].

Conclusions {#S:Conclusion}
===========

We have extended the alternating maximization method (AMM) and modified alternating maximization method (MAMM) given in [@FMPS13] for the computation of best rank one approximation to best ${\mathbf{r}}$-approximations. We have also presented new algorithms such as $2$-alternating maximization method variant (2AMMV) and Newton method for best ${\mathbf{r}}$-approximation (Newton-2). We have provided closed form solutions for computing the $DF$ matrix in Newton-2. We implemented Newton-1 for best rank one approximation [@FMPS11] and Newton-2. From the simulations, we have found out that for rank one approximation of both large and small sized tensors, Newton-1 performed the best. For higher rank approximation, the best performers were 2AMMV in case of large size tensors and AMM or MAMM in case of small size tensors.\
*Acknowledgement*: We thank Daniel Kressner for his remarks.

[MMM]{}

D. Achlioptas and F. McSherry, Fast Computation of Low Rank Approximations, *Proceedings of the 33rd Annual Symposium on Theory of Computing*, 2001, 1–18.

G. Bergqvist, Exact probabilities for typical ranks of $2\times 2\times 2$ and $3\times 3\times 2$ tensors, *Linear Algebra Appl.* 438 (2013), 663–667.

L. de Lathauwer, B. de Moor and J. Vandewalle, A multilinear singular value decomposition, *SIAM J. Matrix Anal. Appl.* 21 (2000), 1253–1278.

L. de Lathauwer, B. de Moor, and J. Vandewalle, On the best rank-1 and rank-$(R_1,R_2, . . . ,R_N)$ approximation of higher-order tensors, *SIAM J. Matrix Anal. Appl.* 21 (2000), pp. 1324–1342.

A. Deshpande and S. Vempala, Adaptive sampling and fast low-rank matrix approximation, Electronic Colloquium on Computational Complexity, Report No. 42 (2006), 1–11.

V. de Silva and L.-H. Lim, Tensor rank and the ill-posedness of the best low-rank approximation problem, *SIAM J. Matrix Anal. Appl.* 30 (2008), 1084–1127.

P. Drineas, R. Kannan and M.W. Mahoney, Fast Monte Carlo algorithms for matrices [I-III]{}: Computing a compressed approximate matrix decomposition, *SIAM J. Comput.* 36 (2006), 132–206.

L. Eldén and B. Savas, A Newton-Grassmann method for computing the best multilinear rank-$(r_1, r_2, r_3)$ approximation of a tensor, *SIAM J. Matrix Anal. Appl.* 31 (2009), 248–271.

S. Friedland, Best rank one approximation of real symmetric tensors can be chosen symmetric, *Front. Math. China* 8 (2013), 19–-40.

S. Friedland, Nonnegative definite hermitian matrices with increasing principal minors, *Special Matrices*, 1 (2013), 1–2.

S. Friedland, *MATRICES*, a book draft in preparation, http://homepages.math.uic.edu/$\sim$friedlan/bookm.pdf, to be published by World Scientific.

S. Friedland, M. Kaveh, A. Niknejad and H. Zare, Fast Monte-Carlo low rank approximations for matrices, *Proc. IEEE Conference SoSE*, Los Angeles, 2006, 218-223.

S. Friedland and V. Mehrmann, Best subspace tensor approximations, arXiv:0805.4220v1.

S. Friedland, V. Mehrmann, A. Miedlar and M. Nkengla, Fast low rank approximations of matrices and tensors, *Journal of Electronic Linear Algebra*, 22 (2011), pp. 1031-1048.

S. Friedland, V. Mehrmann, R. Pajarola and S.K. Suter, On best rank one approximation of tensors, http://arxiv.org/pdf/1112.5914v1.pdf

S. Friedland, V. Mehrmann, R. Pajarola and S.K. Suter, On best rank one approximation of tensors, *Numer. Linear Algebra Appl.* 20 (2013), 942–955.

S. Friedland and G. Ottaviani, The number of singular vector tuples and uniqueness of best rank one approximation of tensors, *Foundations of Computational Mathematics* 14 (2014), 1209–1242.

S. Friedland and A. Torokhti. Generalized rank-constrained matrix approximations, *[SIAM]{} J. Matrix Anal. Appl.* 29 (2007), 656 – 659.

A. Frieze, R. Kannan and S. Vempala, Fast Monte-Carlo algorithms for finding low-rank approximations, *Journal of the ACM*, 51 (2004), 1025–1041.

G.H. Golub and C.F. Van Loan, [*Matrix Computation*]{}, John Hopkins Univ. Press, 3rd Ed., 1996.

S.A. Goreinov, I.V. Oseledets, D.V. Savostyanov, E.E. Tyrtyshnikov and N.L Zamarashkin, How to find a good submatrix, *Matrix methods: theory, algorithms and applications*, 247–256, World Sci. Publ., Hackensack, NJ, 2010.

S.A. Goreinov, E.E. Tyrtyshnikov and N.L. Zamarashkin, A theory of pseudo-skeleton approximations of matrices, *Linear Algebra Appl.* 261 (1997), 1–21.

S.A. Goreinov and E.E. Tyrtyshnikov, The maximum-volume concept in approximation by low-rank matrices, *Contemporary Mathematics* 280 (2001), 47-51.

L. Grasedyck, Hierarchical singular value decomposition of tensors, *SIAM J. Matrix Anal. Appl.* 31 (2010), 2029–2054.

L. Grasedyck, D. Kressner and C. Tobler, A literature survey of low-rank tensor approximation techniques, *GAMM-Mitteilungen* 36 (2013), 53-78.

W. Hackbusch, Tensorisation of vectors and their efficient convolution, *Numer. Math.* 119 (2011), 465–488.

W. Hackbusch, *Tensor Spaces and Numerical Tensor Calculus*, Springer, Heilderberg , 2012.

C.J. Hillar and L.-H. Lim, Most tensor problems are NP-hard, *Journal of the ACM*, 60 (2013), Art. 45, 39 pp.

R.A. Horn and C.R. Johnson, *Matrix Analysis*, Cambridge University Press, 1988.

B.N. Khoromskij, Methods of Tensor Approximation for Multidimensional Operators and Functions with Applications, *Lecture at the workschop*, Schnelle Löser für partielle Differentialgleichungen, Oberwolfach, 18.-23.05, 2008.

T.G. Kolda and B.W. Bader, Tensor decompositions and applications, SIAM Review 51 (2009), pp. 455–500.

R.B. Lehoucq, D.C. Sorensen and C. Yang, *ARPACK User’s Guide : Solution of Large-Scale Eigenvalue Problems With Implicitly Restarted Arnoldi Methods (Software, Environments, Tools)*, SIAM Publications, 1998.

L.-H. Lim, Singular values and eigenvalues of tensors: a variational approach, *Proc. IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing* (CAMSAP ’05), 1 (2005), 129-132.

M.W. Mahoney, M. Maggioni and P. Drineas, Tensor-CUR decompositions for tensor-based data, *Proceedings of the 12th Annual ACM SIGKDD Conference*, 2006, 327–336.

I.V. Oseledets, On a new tensor decomposition, *Dokl. Math.* 80 (2009), 495–496.

I.V. Oseledets, Tensor-Train decompositions, *SIAM J. Sci. Comput.* 33 (2011), 2295–2317.

I.V. Oseledets and E.E. Tyrtyshnikov, Breaking the curse of dimensionality, or how to use SVD in many dimensions, *SIAM J. Sci. Comput.* 31 (2009), 3744–3759.

M. Rudelson and R. Vershynin, Sampling from large matrices: An approach through geometric functional analysis, *Journal of the ACM* 54 (2007), Art. 21, 19 pp.

T. Sarlos, Improved Approximation Algorithms for Large Matrices via Random Projections, Proceedings of the 47th Annual IEEE Symposium on Foundations of Computer Science (FOCS), 2006, 143–152.

B. Savas and L.-H. Lim, Quasi-Newton methods on Grassmannians and multilinear approximations of tensors, *SIAM J. Sci. Comput.* 32 (2010), 3352–3393.

G.W. Stewart On the early history of the singular value decomposition, *SIAM Rev.* 35 (1993), 551–566.

J.M.F. ten Berge and H.A.L. Kiers, Simplicity of core arrays in three-way principal component analysis and the typical rank of $p\times q \times 2$ arrays, *Linear Algebra Appl.* 294 (1999) 169–179.

L.R. Tucker. Some mathematical notes on three-mode factor analysis, *Psychometrika* 31 (1966), 279 – 311.

T. Zhang and G.H. Golub. Rank-one approximation to high order tensors. [*SIAM J. Matrix Anal. Appl*]{}. 23 (2001), pp. 534–550.
