Low-Rank Covariance Function Estimation for Multidimensional Functional Data

First, we give a brief introduction to the popular Tucker decomposition (Tucker, 1966) for finite-dimensional tensors. Let $\mathcal{G}=\bigotimes_{k=1}^{d}{\mathcal{G}}_{k}$ denote a generic tensor product space with finite-dimensional $\mathcal{G}_{k},k=1,\ldots,d$. If the dimension of $\mathcal{G}_{k}$ is $q_{k}$, $k=1,\ldots,d$, then each element in $\mathcal{G}=\bigotimes_{k=1}^{d}{\mathcal{G}}_{k}$ can be identified by an array in $\mathbb{R}^{\prod_{j=k}^{d}q_{k}}$, which contains the coefficients through an orthonormal basis. By Tucker decomposition, any array in $\mathbb{R}^{\prod^{d}_{k=1}q_{k}}$ can be represented in terms of $n$-mode products as follows.

Figure 1 provides a pictorial illustration of a Tucker decomposition. Unlike matrices, the concept of rank is more complicated for arrays of order 3 or above. Tucker decomposition naturally leads to a particular form of rank, called “multilinear rank”, which is directly related to the familiar concept of matrix ranks. To see this, we employ a reshaping operation called matricization, which rearranges elements of an array into a matrix.

For any $\bm{A}\in\mathbb{R}^{q_{1}\times q_{2}\times\cdots\times q_{d}}$, by simple derivations, one can obtain a useful relationship between the $n$-mode matricization and Tucker decomposition $\bm{A}=\bm{G}\times_{1}\bm{U}_{1}\times_{2}\cdots\times_{d}\bm{U}_{d}$:

$$\bm{A}_{(n)}=\bm{U}_{n}\bm{G}_{(n)}(\bm{U}_{d}\otimes\cdots\otimes\bm{U}_{n+1}\otimes\bm{U}_{n-1}\otimes\cdots\otimes\bm{U}_{1})^{\intercal},$$

(4)

where, with a slight abuse of notation, $\otimes$ also represents the Kronecker product between matrices. Hence if the factor matrices are of full rank, then $\mathrm{rank}(\bm{A}_{(n)})=\mathrm{rank}(\bm{G}_{(n)})$. The vector $(\mathrm{rank}(\bm{A}_{(1)}),\dots,\mathrm{rank}(\bm{A}_{(d)}))$ is known as the multilinear rank of $\bm{A}$. Clearly from (4), one can choose a Tucker decomposition such that $\{\bm{U}_{k}:k=1,\ldots,d\}$ are orthonormal matrices and $\mathrm{rank}(\bm{U}_{k})=r_{k}$. Therefore a “small” multilinear rank corresponds to a small core tensor and thus an intrinsic dimension reduction, which potentially improves estimation and interpretation. We will relate this low-rank structure to multidimensional functional data.

To encourage low-rank structures in covariance function estimation, we generalize the matricization operation for finite-dimensional arrays to infinite-dimensional tensors (Hackbusch, 2012). Here let $\mathcal{G}=\bigotimes_{k=1}^{d}{\mathcal{G}}_{k}$ denote a generic tensor product space where $\mathcal{G}_{k}$ is an RKHS of functions with an inner product $\langle\cdot,\cdot\rangle_{{\mathcal{G}}_{k}}$, for $k=1,\ldots,d$.

Notice that the tensor product space ${\mathcal{G}}=\bigotimes_{k=1}^{d}{\mathcal{G}}_{k}$ can be generated by some elementary tensors of the form $\bigotimes^{d}_{k=1}f_{k}(x_{1},\dots,x_{d})=\prod^{d}_{k=1}f_{k}(x_{k})$ where $f_{k}\in{\mathcal{G}}_{k},k=1,\ldots,d$. More specifically, $\mathcal{G}$ is the completion of the linear span of all elementary tensors under the inner product $\langle\bigotimes^{d}_{k=1}f_{k},\bigotimes^{d}_{k=1}f^{\prime}_{k}\rangle_{\mathcal{G}}=\prod^{d}_{k=1}\langle f_{k},f_{k}^{\prime}\rangle_{{\mathcal{G}}_{k}}$, for any $f_{k},f_{k}^{\prime}\in{\mathcal{G}}_{k}$.

In Definition 4 below, we generalize matricization/unfolding for finite-dimensional arrays to infinite-dimensional elementary tensors. We also define a square unfolding for infinite-dimensional tensors that will be used to describe the spectrum of covariance operators.

Note that the range of each functional unfolding operator, either ${\mathcal{U}}_{j}$, $j=1,\ldots,d$ or $\mathcal{S}$, is a tensor product of two RKHSs, so its output can be interpreted as an (induced) operator. Given a function $f\in\mathcal{G}$, the multilinear rank can be defined as $(\mathrm{rank}(f_{(1)}),\dots,\mathrm{rank}(f_{(d)}))$, where $f_{(j)}$’s are interpreted as an operator here and $\mathrm{rank}(A)$ is the rank of any operator $A$. If all $\mathcal{G}_{k},k=1,\ldots,d$ are finite-dimensional, the singular values of the output of any functional unfolding operator match with those of the $j$-mode matricization (of the corresponding array representation).

Suppose that the random field $X\in{\mathcal{H}}=\bigotimes_{k=1}^{p}{\mathcal{H}}_{k}$ where each ${\mathcal{H}}_{k}$ is a RKHS of functions equipped with an inner product $\langle\cdot,*\rangle_{k}$ and corresponding norm $\|\cdot\|_{k}$, $k=1,\ldots,p$. Then the covariance function $\gamma_{0}$ resides in ${\mathcal{H}}\otimes{\mathcal{H}}=(\bigotimes_{j=1}^{p}{\mathcal{H}}_{j})\otimes(\bigotimes_{k=1}^{p}{\mathcal{H}}_{k})$. To estimate $\gamma_{0}$, we could consider a special case of $\mathcal{G}=\bigotimes^{d}_{j=1}\mathcal{G}_{j}$ in Section 3.2 by letting $d=2p$, $\mathcal{G}_{j}=\mathcal{H}_{j}$ for $j=1,\dots,p$; $\mathcal{G}_{j}=\mathcal{H}_{j-p}$ for $j=p+1,\dots,d$; and $\langle\cdot,\cdot\rangle_{{\mathcal{G}}_{j}}=\langle\cdot,\cdot\rangle_{j}$ for $j=1,\ldots,d$.

Clearly, the elements of ${\mathcal{H}}\otimes{\mathcal{H}}$ are identified by those in $\mathcal{G}=\bigotimes^{d}_{j=1}\mathcal{G}_{j}$. In terms of the folding structure, ${\mathcal{H}}\otimes{\mathcal{H}}$ has a squarely unfolded structure. Since a low-multilinear-rank structure is represented by different unfolded forms, it would be easier to study the completely folded space $\bigotimes^{d}_{k=1}{\mathcal{G}}_{k}$ instead of the squarely unfolded space ${\mathcal{H}}\otimes{\mathcal{H}}$. We use $\Gamma_{0}$ to represent the folded covariance function, the corresponding element of $\gamma_{0}$ in ${\mathcal{G}}$. In other words, $\Gamma_{0,\mathrel{\scalebox{0.5}{$\blacksquare$}}}=\gamma_{0}$. For any $\Gamma\in{\mathcal{G}}$, $\mathrm{rank}(\Gamma_{\mathrel{\scalebox{0.5}{$\blacksquare$}}})$ is defined as the two-way rank of $\Gamma$ while $\mathrm{rank}(\Gamma_{(1)}),\dots,\mathrm{rank}(\Gamma_{(p)})$ are defined as the one-way ranks of $\Gamma$.

Here we illustrate the roles of one-way and two-way ranks in the modeling of covariance functions. For a general $\mathcal{G}=\bigotimes^{d}_{j=1}\mathcal{G}_{j}$, let $\{e_{k,l_{k}}:l_{j}=1,\dots,q_{k}\}$ be a set of orthonormal basis functions of ${\mathcal{G}}_{k}$ for $k=1,\dots,d=2p$, where $q_{k}$ is allowed to be infinite, depending on the dimensionality of ${\mathcal{G}}_{k}$. Then $\{\bigotimes^{d}_{k=1}e_{k,l_{k}}:l_{k}=1,\dots,q_{k};k=1,\dots,d\}$ forms a set of orthonormal basis functions for ${\mathcal{G}}$. Thus for any $\Gamma\in{\mathcal{G}}$, we can express

$\displaystyle\Gamma$

$\displaystyle=\sum_{k_{1},k_{2},\dots,k_{d}}\bm{B}_{k_{1},\dots,k_{d}}\bigotimes_{i=1}^{d}e_{i,k_{i}},$

(5)

where the coefficients $\bm{B}_{k_{1},\dots,k_{d}}$ are real numbers. For convenience, we collectively put them into an array $\bm{B}\in\mathbb{R}^{\prod^{d}_{k=1}q_{k}}$.

To illustrate the low-multilinear-rank structures for covariance functions, we consider $p=2$, i.e., $d=2p=4$, and then by (5) the folded covariance function $\Gamma$ can be expressed by

$$\Gamma(s_{1},s_{2},t_{1},t_{2})=\sum_{k_{1}=1}^{q_{1}}\sum_{k_{2}=1}^{q_{2}}\sum_{k_{3}=1}^{q_{1}}\sum_{k_{4}=1}^{q_{2}}\bm{B}_{k_{1},k_{2},k_{3},k_{4}}e_{1,k_{1}}(s_{1})e_{2,k_{2}}(s_{2})e_{1,k_{3}}(t_{1})e_{2,k_{4}}(t_{2}).$$

To be precise, the covariance function is the squarely unfolded $\Gamma_{{\mathrel{\scalebox{0.5}{$\blacksquare$}}}}((s_{1},s_{2}),(t_{1},t_{2}))\equiv\Gamma(s_{1},s_{2},t_{1},t_{2})$. Suppose that $\bm{B}$ possesses (or is well-approximated by) a structure of a low multilinear rank, and yields Tucker decomposition $\bm{B}=\bm{E}\times_{1}\bm{U}_{1}\times_{2}\bm{U}_{2}\times_{3}\bm{U}_{1}\times_{4}\bm{U}_{2}$²²2Definition 1 is extended to the case when $q_{n}$ is infinite. where $\bm{E}\in\mathbb{R}^{r_{1}\times r_{2}\times r_{1}\times r_{2}}$, $\bm{U}_{k}\in\mathbb{R}^{q_{k}\times r_{k}}$ for $k=1,2$, and columns of $\bm{U}_{k}$ are orthonormal. Apparently $R:=\mathrm{rank}(\bm{B}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}})$ is the two-way rank of $\Gamma$, while $r_{1}$ and $r_{2}$ are the corresponding one-way ranks. Now the covariance function can be further represented as

$$\Gamma(s_{1},s_{2},t_{1},t_{2})=\sum_{j_{1}=1}^{r_{1}}\sum_{j_{2}=1}^{r_{2}}\sum_{j_{3}=1}^{r_{1}}\sum_{j_{4}=1}^{r_{2}}\bm{E}_{j_{1},j_{2},j_{3},j_{4}}u_{j_{1}}(s_{1})v_{j_{2}}(s_{2})u_{j_{3}}(t_{1})v_{j_{4}}(t_{2}),$$

where $\{u_{j}:j=1,\ldots,r_{1}\}$ and $\{v_{k}:k=1,\ldots,r_{2}\}$ are (possibly infinite) linear combinations of the original basis functions. In fact, $\{u_{j}:j=1,\ldots,r_{1}\}$ and $\{v_{k}:k=1,\ldots,r_{2}\}$ are the sets of orthonormal functions of ${\mathcal{G}}_{1}$ and ${\mathcal{G}}_{2}$ respectively. Apparently $\mathrm{rank}(\bm{E}_{{\mathrel{\scalebox{0.5}{$\blacksquare$}}}})=R$.

Consider the eigen-decomposition of the squarely unfolded core tensor $\bm{E}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}=\bm{P}\bm{D}\bm{P}^{T}$ where $\bm{D}=\mathrm{diag}(\lambda_{1},\lambda_{2},\dots,\lambda_{R})$ and $\bm{P}\in\mathbb{R}^{r_{1}r_{2}\times R}$ has orthonormal columns. Then we obtain the eigen-decomposition of the covariance function $\Gamma_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}$:

$$\Gamma_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}((s_{1},s_{2}),(t_{1},t_{2}))=\sum_{g=1}^{R}\lambda_{g}f_{g}(s_{1},s_{2})f_{g}(t_{1},t_{2}),$$

where the eigenfunction is

$$f_{g}(s_{1},s_{2})=\sum_{j_{1}=1}^{r_{1}}\sum_{j_{2}=1}^{r_{2}}\bm{P}_{j_{2}+(j_{1}-1)r_{1},g}u_{j_{1}}(s_{1})v_{j_{2}}(s_{2})=:\,\begin{cases}\sum^{r_{1}}_{j_{1}=1}a_{j_{1},g}(s_{2})u_{j_{1}}(s_{1})\\ \sum^{r_{2}}_{j_{2}=1}b_{j_{2},g}(s_{1})v_{j_{2}}(s_{2})\end{cases},$$

with $a_{j_{1},g}(\cdot)=\sum_{j_{2}=1}^{r_{2}}\bm{P}_{j_{2}+(j_{1}-1)r_{1},g}v_{j_{2}}(\cdot)$ and $b_{j_{2},g}(\cdot)=\sum_{j_{1}=1}^{r_{1}}\bm{P}_{j_{2}+(j_{1}-1)r_{1},g}u_{j_{1}}(\cdot)$.

First, this indicates that the two-way rank $R$ is the same as the rank of the covariance operator. Second, this shows that $\{u_{j_{1}}:j_{1}=1,\ldots,r_{1}\}$ is the common basis for the variation along the dimension $s_{1}$, hence describing the marginal structure along $s_{1}$. Similarly $\{v_{j_{2}}:j_{2}=1,\ldots,r_{2}\}$ is the common basis that characterizes the marginal variation along the dimension $s_{2}$. We call them the marginal basis along the respective dimension. Therefore, the one-way ranks $r_{1}$ and $r_{2}$ are the minimal numbers of the one-dimensional functions for the dimensions $s_{1}$ and $s_{2}$ respectively that construct all the eigenfunctions of covariance function $\Gamma$.

Similarly, for $p$-dimensional functional data, each eigenfunction can be represented by a linear combination of $p$-products of univariate functions. One can then show that the two-way rank $R$ is the same as the rank of the covariance operator and the one-way ranks $r_{1},\dots,r_{p}$ are the minimal numbers of one-dimensional functions along respective dimensions that characterize all eigenfunctions of the covariance operator.

Compared to typical low-rank covariance modelings only in terms of $R$, we also intend to regularize the one-way ranks $r_{1},\dots,r_{p}$ for two reasons. First, the illustration above shows that the structure of low one-way ranks encourages a “sharing” structure of one-dimensional variations among different eigenfunctions. Promoting low one-way ranks can facilitate additional dimension reduction and further alleviates the curse of dimensionality. Moreover, one-dimensional marginal structures will provide more details of the covariance function structure and thus help with a better understanding of $p$-dimensional eigenfunctions.

Therefore, we will utilize both one-way and two-way structures and propose an estimation procedure that regularizes one-way and two-way ranks jointly and flexibly, with the aim of seeking the “sharing” of marginal structures while controlling the number of eigen-components simultaneously.

Before deriving a computational algorithm, we notice that the optimization (7) is an infinite-dimensional optimization which is generally unsolvable. To overcome this challenge, we show that the solution to (7) always lies in a known finite-dimensional sub-space given data, hence allowing a finite-dimensional parametrization. Indeed, we are able to achieve a stronger result in Theorem 1 which holds for the general class of estimators obtained by (6).

Let $\mathcal{L}_{n,m}=\left\{T_{ijk}:i=1,\dots,n,j=1,\dots,m,k=1,\dots,p\right\}$.

The proof of Theorem 1 is given in Section S1 of the supplementary material. By Theorem 1, we can now only focus on covariance function estimators of the form (9). Let $\bm{B}=\bm{A}\times_{1}\bm{M}_{1}^{T}\cdots\times_{p}\bm{M}_{p}^{T}\times_{p+1}\bm{M}_{1}^{T}\cdots\times_{2p}\bm{M}_{p}^{T}$, where $\bm{M}_{k}$ is a $nm\times q_{k}$ matrix such that $\bm{M}_{k}\bm{M}_{k}^{T}=\bm{K}_{k}=\left[K(T_{i_{1},j_{1},k},T_{i_{2},j_{2},k})\right]_{1\leq i_{1},i_{2}\leq n,1\leq j_{1},j_{2}\leq m}$. With $\bm{B}$, we can express

$\displaystyle\begin{split}\Gamma(s_{1},\ldots,s_{p},t_{1},\ldots,t_{p})=\bm{B}&\times_{1}\{\bm{M}_{1}^{+}\bm{z}_{1}(s_{1})\}^{\intercal}\cdots\times_{p}\{\bm{M}_{p}^{+}\bm{z}_{p}(s_{p})\}^{\intercal}\\ &\times_{p+1}\{\bm{M}_{1}^{+}\bm{z}_{1}(t_{1})\}^{\intercal}\cdots\times_{2p}\{\bm{M}_{p}^{+}\bm{z}_{p}(t_{p})\}^{\intercal},\end{split}$

(10)

where $z_{k}(\cdot)$ is defined in Theorem 1 and $\bm{M}_{k}^{+}$ is the Moore–Penrose inverse of matrix $\bm{M}_{k}$.

The Gram matrix $\bm{K}_{k}$ is often approximately low-rank. For computational simplicity, one could adopt $q_{k}$ to be significantly smaller than $nm$. Ideally we can obtain the “best” low-rank approximation with respect to the Frobenius norm by eigen-decomposition, but a full eigen-decomposition is computationally expensive. Instead, randomized algorithms can be used to obtain low-rank approximations in an efficient manner (Halko et al., 2009).

One can easily show that the eigenvalues of the operator $\Gamma_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}$ are the same as those of the matrix ${\bm{B}}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}$ and that the singular values of the operator $\Gamma_{(j)}$ are the same as those of the matrix ${\bm{B}}_{(j)}$. Therefore, solving (7) is equivalent to solving the following optimization:

$$\min_{\bm{B}}\left\{\tilde{\ell}_{\mathrm{square}}(\bm{B})+\lambda\left[\beta h({\bm{B}}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}})+\frac{1-\beta}{p}\sum_{k=1}^{p}\left\|{\bm{B}}_{(j)}\right\|_{*}\right]\right\},$$

(11)

where $\|\cdot\|_{*}$ also represents the trace norm of matrices, $h(\bm{H})=\left\|\bm{H}\right\|_{*}$ if matrix $\bm{H}$ is positive semi-definite, and $h(\bm{H})=\infty$ otherwise, and $\tilde{\ell}_{\mathrm{square}}(\bm{B})=\ell_{\mathrm{square}}(\Gamma)$, where $\Gamma$ is constructed from (10).

Beyond estimating the covariance function, one may be further interested in the eigen-decomposition of $\Gamma_{{\mathrel{\scalebox{0.5}{$\blacksquare$}}}}$ via the $L^{2}$ inner product, e.g., to perform functional principal component analysis in the usual sense. Due to the finite-dimensional parametrization, a closed-form expression of $L^{2}$ eigen-decomposition can be derived from our estimator without further discretization or approximation. In addition, we can obtain a similar one-way analysis in terms of the $L_{2}$ inner product. We can define a $L^{2}$ singular value decomposition via the Tucker form and obtain the $L^{2}$ marginal basis. Details are given in Appendix A.

We solve (11) by the accelerated alternating direction method of multipliers (ADMM) algorithm (Kadkhodaie et al., 2015). We begin with an alternative form of (11):

		$\displaystyle\min_{\bm{B}\in\mathbb{R}^{q_{1}\times\cdots\times q_{2p}}}\left\{\tilde{\ell}_{\mathrm{square}}(\bm{B})+\lambda\beta h(\bm{D}_{0,{\mathrel{\scalebox{0.5}{$\blacksquare$}}}})+\lambda\frac{1-\beta}{p}\sum_{k=1}^{p}\left\|\bm{D}_{j,(j)}\right\|_{*}\right\}.$		(12)
		$\displaystyle\text{subject to}\ \ \ \bm{B}=\bm{D}_{0}=\bm{D}_{1}=\cdots=\bm{D}_{p}$		(13)

where $q_{p+k}=q_{k}$ for $k=1,\dots,p$.

Then a standard ADMM algorithm solves the optimization problem (12) by minimizing the augmented Lagrangian with respect to different variables alternatively. More explicitly, at the $(t+1)$-th iteration, the following updates are implemented:

	$\displaystyle\bm{B}^{(t+1)}$	$\displaystyle=\operatorname*{argmin}_{\bm{B}}\left\{\tilde{\ell}_{\mathrm{square}}(\bm{B})+\frac{\eta}{2}\|\bm{B}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}-\bm{D}_{0,{\mathrel{\scalebox{0.5}{$\blacksquare$}}}}^{(t)}+\bm{V}_{0,{\mathrel{\scalebox{0.5}{$\blacksquare$}}}}^{(t)}\|_{F}^{2}+\frac{\eta}{2}\sum_{k=1}^{p}\left\|{\bm{B}}_{(k)}-\bm{D}_{k,(k)}^{(t)}+\bm{V}_{k,(k)}^{(t)}\right\|_{F}^{2}\right\},$		(14a)
	$\displaystyle\bm{D}_{0}^{(t+1)}$	$\displaystyle=\operatorname*{argmin}_{\bm{D}_{0}}\left\{\lambda\beta h(\bm{D}_{0,{\mathrel{\scalebox{0.5}{$\blacksquare$}}}})+\frac{\eta}{2}\left\|\bm{B}^{(t+1)}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}-\bm{D}_{0,{\mathrel{\scalebox{0.5}{$\blacksquare$}}}}+\bm{V}_{0,{\mathrel{\scalebox{0.5}{$\blacksquare$}}}}^{(t)}\right\|_{F}^{2}\right\},$		(14b)
	$\displaystyle\bm{D}_{k}^{(t+1)}$	$\displaystyle=\operatorname*{argmin}_{\bm{D}_{k}}\left\{\lambda\frac{1-\beta}{p}\|\bm{D}_{k,(k)}\|_{*}+\frac{\eta}{2}\left\|\bm{B}^{(t+1)}_{(k)}-\bm{D}_{k,(k)}+\bm{V}_{k,(k)}^{(t)}\right\|_{F}^{2}\right\},\ k=1,\dots,p,$		(14c)
	$\displaystyle\bm{V}_{k}^{(t+1)}$	$\displaystyle=\bm{V}_{k}^{(t)}+\bm{B}^{(t+1)}-\bm{D}_{k}^{(t+1)},\ k=0,\dots,p,$		(14d)

where $\bm{V}_{k}\in\mathbb{R}^{q_{1}\times\cdots q_{2p}}$, for $k=0,\dots,p$, are scaled Lagrangian multipliers and $\eta>0$ is an algorithmic parameter. An adaptive strategy to tune $\eta$ is provided in Boyd et al. (2010). One can see that Steps (14a), (14b) and (14c) involve additional optimizations. Now we discuss how to solve them.

The objective function of (14a) is a quadratic function, and so we can easily solve this with a closed-form solution, given in line 2 of Algorithm 1. To solve (14b) and (14c), we use proximal operator $\mathrm{prox}^{k}_{v}$, $k=1,\dots,p$ and $\mathrm{prox}_{v}^{+}:\mathbb{R}^{q_{1}\times\cdots\times q_{2p}}\rightarrow\mathbb{R}^{q_{1}\times\cdots\times q_{2p}}$ respectively defined by

	$\displaystyle\mathrm{prox}_{v}^{k}(\bm{A})$	$\displaystyle=\operatorname*{argmin}_{\bm{W}\in\mathbb{R}^{q_{1}\times\cdots\times q_{2p}}}\left\{\frac{1}{2}\|\bm{W}_{(k)}-\bm{A}_{(k)}\|_{F}^{2}+v\|\bm{W}_{(k)}\|_{*}\right\},$		(15a)
	$\displaystyle\mathrm{prox}^{+}_{v}(\bm{A})$	$\displaystyle=\operatorname*{argmin}_{\bm{W}\in\mathbb{R}^{q_{1}\times\cdots\times q_{2p}}}\left\{\frac{1}{2}\|\bm{W}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}-\bm{A}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}\|_{F}^{2}+vh(\bm{W}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}})\right\},$		(15b)

for $v\geq 0$. By Lemma 1 in Mazumder et al. (2010), the solutions to (15) have closed forms.

For (15a), write the singular value decomposition of $\bm{A}_{(k)}$ as $\bm{U}\mathrm{diag}((\tilde{a}_{1},\dots,\tilde{a}_{q_{k}}))\bm{V}^{\intercal}$, then $[\mathrm{prox}_{v}^{k}(\bm{A})]_{(k)}=\bm{U}\mathrm{diag}(\bm{\tilde{}}{\bm{c}})\bm{V}^{\intercal}$ where $\tilde{\bm{c}}=((\tilde{a}_{1}-v)_{+},(\tilde{a}_{2}-v)_{+},\dots,(\tilde{a}_{q_{k}}-v)_{+})$. As for (15b), is restricted to be a symmetric matrix since the penalty $h$ equals infinity otherwise. Thus (15b) is equivalent to minimizing $\left\{(1/2)\|\bm{W}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}-(\bm{A}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}+\bm{A}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}^{\intercal})/2\|_{F}^{2}+vh(\bm{W}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}})\right\}$ since $\langle\bm{W}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}},(\bm{A}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}-\bm{A}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}^{\intercal})/2\rangle=\langle(\bm{W}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}+\bm{W}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}^{\intercal})/2,(\bm{A}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}-\bm{A}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}^{\intercal})/2\rangle=0$ . Suppose that $(\bm{A}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}+\bm{A}_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}^{\intercal})/2$ yields eigen-decomposition $\bm{P}\mathrm{diag}((\tilde{a}_{1},\dots,\tilde{a}_{q})\bm{P}^{\intercal}$. Then $[\mathrm{prox}^{+}_{v}(\bm{A})]_{\mathrel{\scalebox{0.5}{$\blacksquare$}}}=\bm{P}\mathrm{diag}(\bm{\tilde{}}{\bm{c}})\bm{P}^{\intercal}$, where $\tilde{\bm{c}}=((\tilde{a}_{1}-v)_{+},(\tilde{a}_{2}-v)_{+},\dots,(\tilde{a}_{q}-v)_{+})$. Unlike singular values, the eigenvalues may be negative. Hence, as opposed to $\mathrm{prox}^{k}_{v}$, this procedure $\mathrm{prox}_{v}^{+}$ also removes eigen-components with negative eigenvalues.