Higher-Order Low-Rank Regression
Guillaume Rabusseau\footnote{Contact author: guillaume.rabusseau@lif.univ-mrs.fr
The purpose of this document is to provide both the basic paper template and submission guidelines. Abstracts should be a single paragraph, between 4--6 sentences long, ideally. Gross violations will trigger corrections at the camera-ready phase.
abstract.tex

This paper proposes an efficient algorithm (HOLRR) to handle regression tasks where
the outputs have a tensor structure. We formulate the regression problem as the minimization 
of a least square criterion under a multilinear rank constraint, a difficult  non convex problem. 
HOLRR computes efficiently an approximate solution of this problem, with solid theoretical guarantees.
A kernel extension is also presented. Experiments on synthetic and real
data show that HOLRR outperforms multivariate and multilinear regression methods and
is considerably faster than existing tensor methods.


%This algorithm tackles the minimization of a 
%least square criterion under a mutlilinear rank constraint. Although this 
%is a difficult non-convex optimization problem, we show that HOLRR provides
%good theoretical approximation guarantees and is computationally efficient. A
\comment{
The predictive accuracy of HOLRR is assessed
on experiments 

Extending univariate and multivariate methods to tensor-structured data is a challenging task which has recently 
received a growing interest.
In this paper, we generalize the reduced-rank regression method from vector to tensor-structured response variables. 
Reduced-rank regression allows one to capture dependencies in the response variables by enforcing the regression
matrix to be of low rank, but it is unable to capture higher-order dependencies. 
To address this limitation we propose to learn a regression tensor of low multilinear
rank. We design an efficient algorithm for which we provide theoretical approximation guarantees and a kernelized extension.
We compare our method  with other multi-output regression methods on both synthetic and real data 
where it exhibits better predictive accuracy than standard multivariate methods and is more computationally
efficient than existing tensor regression methods.}


 
algo.tex

\section{Low-Rank Regression for Tensor-Valued Functions}
\label{sec:algo}
\subsection{Problem Formulation}


We consider a multivariate regression task where the response has a tensor structure. 
Let  $f: \R^{d_0}\to\R^{d_1\times d_2\times\cdots\times d_p}$ be the function
we want to learn from a sample of input-output data $\{(\vec{x}\nsample{n}, \ten{Y}\nsample{n})\}_{n=1}^N$ 
drawn from the model  $\ten{Y} = f(\vec{x}) + \ten{E}$, where $\ten{E}$ is an error term.
We assume that the function $f$ is linear,
that is $f(\vec{x}) = \ten{W} \ttv{1} \vec{x}$ for some regression tensor 
$\ten{W} \in  \R^{d_0\times d_1 \times \cdots \times d_p}$. 
Note that the vectorization of this relation leads to $\vectorize{f(\vec{x})} = \tenmat{W}{1}^\top \vec{x}$
showing that this model is equivalent to the standard multivariate linear model.

%
\paragraph{Vectorization of the outputs.}
One way to tackle this linear regression task for tensor responses would be to vectorize 
each output sample and  to perform a standard multivariate low-rank regression 
on the data $\{(\vec{x}\nsample{n}, \vectorize{\Y\nsample{n}})\}_{n=1}^N\subset \R^{d_0}\times\R^{d_1\cdots d_p}$. 
%
%A major drawback of this approach is that the tensor structure of the output is completely lost because of the vectorization step. It can capture linear dependencies between components of the output tensor but neglects the coupling between the modes of the tensor: . For illustration, suppose that the output matrix encodes the samples of $d_1$ continuous variables at $d_2$ different time steps.  Concatenating the samples of all the $d_1$ variables in a vector of size $d_1\times d_2$ and applying multivariate regression cannot, for instance, make any distinction between output samples that belong to different variables with those that are inside the same variable, so it is not surprising that standard multivariate~(vector-valued) regression analysis cannot model high level output structure that one could expect from a tensor~(matrix) representation of the output.
A major drawback of this approach is that the tensor structure of the output is lost in the vectorization step. 
The low-rank model tries to capture linear dependencies between components of the output 
but it ignores \emph{higher level dependencies} that could be present in a tensor-structured output. 
For illustration, suppose the output is a matrix encoding the samples of $d_1$ continuous variables at $d_2$ different time steps,
 one could expect structural relations between the $d_1$ time series, for example linear dependencies between the rows 
 of the output matrix.
  

\paragraph{Low-rank regression for tensor responses.} 
To overcome the limitation described above we propose an  extension of the low-rank regression method 
 for tensor-structured responses by enforcing low multilinear rank of the regression tensor $\ten{W}$. 
Let $\{(\vec{x}\nsample{n}, \ten{Y}\nsample{n})\}_{n=1}^N 
\subset \R^{d_0}\times \R^{d_1\times d_2\times \cdots\times d_p}$ be  a  training sample of input/output data
drawn from the model $f(\vec{x}) = \ten{W} \ttv {1}\vec{x} + \ten{E}$ where $\ten{W}$ is assumed of low
multilinear rank. 
Considering the framework of empirical risk minimization, we want to find a low-rank regression tensor 
$\W$ minimizing the loss on the training data. To avoid numerical 
instabilities and to prevent overfitting we add a ridge regularization to the objective
function, leading to the following minimization problem
\begin{align}
\label{eq:holrr-general-loss}
\MoveEqLeft\min_{\W\in \R^{d_0\times  \cdots \times d_p}} \sum_{n=1}^N
\ell(\W\ttv{1}\vec{x}\nsample{n}, \ten{Y}\nsample{n}) + \gamma \norm{\ten{W}}^2_F \\
 &\text{ s.t. } \rank(\W) \leq (R_0,R_1,\cdots,R_p)\nonumber\enspace,
\end{align}
for some given integers $R_0,R_1,\cdots,R_p$ and where $\ell$ is a loss function. 
In this paper, we consider the squared error loss
between tensors defined by
$\ell(\ten{T},\hat{\ten{T}}) = \norm{\ten{T}-\hat{\ten{T}}}^2_F$.  
%
Using this loss we can rewrite problem~(\ref{eq:holrr-general-loss}) as
\begin{align}\label{eq:holrr-pbm}
\MoveEqLeft\min_{\W\in \R^{d_0\times d_1 \times \cdots \times d_p}} \norm{\W\ttm{1}\X - \Y}^2_F 
+ \gamma \norm{\W}^2_F \\
&\text{s.t. }\rank(\W) \leq  (R_0, R_1, \cdots, R_p)\enspace,\nonumber
\end{align}
where the input matrix $\X\in \R^{N\times d_0}$ and the output tensor
$\Y\in \R^{N\times d_1 \times \cdots \times d_p}$ are defined by
$\X_{n,:} = (\vec{x}\nsample{n})^\top$, $\Y_{n,:,\cdots,:} = \Y\nsample{n}$
for $n=1,\cdots, N$ ($\Y$ is the tensor obtained by stacking the output
tensors along the first mode).

\paragraph{Low-rank regression function.} Let $\W^*$ be a solution of problem~(\ref{eq:holrr-pbm}), it
follows from the multilinear rank constraint that $\W^* = \G \ttm{1}\U_0\ttm{2}\cdots\ttm{p+1}\U_p$
for some core tensor $\G \in \R^{R_0\times\cdots\times R_p}$ and orthogonal matrices 
$\U_i\in\R^{d_i\times R_i}$ for $0\leq i \leq p$. The regression function $f^* : \vec{x} \mapsto \W^*\ttv{1}\vec{x}$
can thus be written as 
$$f^*:\x \mapsto \G \ttm{1}\x^\top\U_0\ttm{2}\cdots\ttm{p+1}\U_p\enspace .$$

This implies several interesting properties. First, for any $\vec{x}\in \R^{d_0}$ we
have $f^*(\x) = \T_\x \ttm{1}\U_1\ttm{2}\cdots\ttm{p} \U_p$ with 
$\T_\x = \G \ttv{1} \U_0^\top\x$, which implies $\rank(f^*(\vec{x})) \leq (R_1,\cdots,R_p)$,
that is the image of $f^*$ is a set of tensors with low multilinear rank. Second, the relation
between $\x$ and $\Y=f^*(\x)$ is explained in a low dimensional subspace of size 
$R_0\times R_1\times\cdots\times R_p$. Indeed one can decompose the mapping $f^*$ into the 
following steps: (i) project $\x$ in $\R^{R_0} $ as $\bar{\x} = \U_0^\top\x$, (ii) perform
a low-dimensional mapping  $\bar{\Y} = \G\ttv{1}\bar{\x}$, (iii) project back into the
output space to get $\Y = \bar{\Y}\ttm{1}\U_1\ttm{2}\cdots\ttm{p}\U_p$.


\paragraph{Comparison with LRR.} First, it is obvious that the function learned by vectorizing
the outputs and performing LRR does not enforce any multilinear structure in the predictions.
That is, the tensor obtained by reshaping a vectorized prediction will not necessarily have a 
low multilinear rank. This is well illustrated in the experiment from Section~\ref{sec:xp_images}
where this fact is particularly striking when we try to learn an image from noisy measurements 
with a rank constraint set to $1$ (see Figure~\ref{fig:xp_images}, left).
Second, since the multilinear rank is defined with respect to the ranks of the matricizations 
there are some scenarios where the rank constraint in LRR after vectorizing the output will 
be (almost) equivalent to the multilinear rank constraint in~(\ref{eq:holrr-pbm}). 
This is clear if the regression tensor is of low mode-$1$ rank and of full mode-$n$ rank for
$n > 1$, but even when $R_0\ll R_i$ for $i\geq 1$ we can expect  LRR to be
able to capture most of the low-rank structure which is concentrated along the first mode.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



\subsection{Higher-Order Low-Rank Regression}
\label{sec:holrr}
We now propose an efficient algorithm to tackle problem~(\ref{eq:holrr-pbm}).
Theoretical approximation guarantees are given in the next section.

%The rank constraint in~(\ref{eq:holrr-pbm}) is equivalent to enforcing the decomposition 
%\begin{equation}
%\W = \G\ttm{1}\U_0\ttm{2}\U_1\cdots\ttm{3}\U_p
%\end{equation}
%for some $\G\in \R^{R_0\times \cdots \times R_p}$ and orthogonal matrices 
%$\U_i\in\R^{d_i\times R_i}$ for $0\leq  i\leq p$. 
\paragraph{Problem reformulation.}  We first show that the ridge regularization term in
~(\ref{eq:holrr-pbm}) can be incorporated in the data fitting term.
Let $\tX\in \R^{(N+d_0)\times d_0}$ and $\tY\in \R^{(N+d_0)\times d_1 \times \cdots \times d_p}$
be defined by $\tX^\top = \left(\X\ |\ \gamma\I\right)^\top$ and 
$\tenmat{\tilde{Y}}{1}^\top = \left(\tenmat{Y}{1}\ |\ \mats{0}\right)^\top$.
It is easy to check that the objective function in~(\ref{eq:holrr-pbm}) is equal to 
$\norm{\W\ttm{1}\tX - \tY}^2_F$. Minimization problem~(\ref{eq:holrr-pbm})
can then be rewritten as
\begin{align}\label{eq:holrr-pbm2}
\MoveEqLeft\hspace{-1.5cm}\min_{\G\in \R^{R_0\times R_1 \times \cdots \times R_p},
\atop \U_i \in \R^{d_i\times R_i} \text{ for } 0\leq i \leq p} 
\norm{\W\ttm{1}\tX - \tY}^2_F  \\
\quad\text{s.t. }&\W = \G\ttm{1}\U_0\ttm{2}\U_1\cdots\ttm{p+1}\U_p,\nonumber\\
\quad\ \ \ \ \ \ &\U_i^\top\U_i=\I \text{  for } 0\leq i \leq p\enspace.\nonumber
\end{align}

We now show that this minimization problem can be reduced to finding $p+1$
projection matrices onto subspaces of dimension $R_0,R_1,\cdots, R_p$.
We start by showing that the core tensor $\G$ solution  of~(\ref{eq:holrr-pbm2})
is determined by the factor matrices $\U_0,\cdots,\U_p$. 
\begin{theorem}\label{thm:core-tensor}
For given orthogonal matrices $\U_0,\cdots,\U_p$ the tensor $\G$ that
minimizes~(\ref{eq:holrr-pbm2}) is given by
$$\G = \tY \ttm{1}(\U_0^\top\tX^\top\tX\U_0)^{-1}\U_0^\top\tX^\top
\ttm{2}\U_1^\top\ttm{3}\cdots\ttm{p+1}\U_p^\top\enspace.$$
%\begin{align*}
%\G = \Y &\ttm{1}(\U_0^\top(\X^\top\X+\gamma\I)\U_0)^{-1}\U_0^\top\X^T\\
%&\ttm{2}\U_1^\top\ttm{3}\cdots\ttm{p+1}\U_p^\top\enspace.
%\end{align*}
\end{theorem}
\begin{proof}
Using Eq.~\ref{eq:vecTucker} the objective function in~(\ref{eq:holrr-pbm2}) can be written as
$$\norm{(\U_p\kron \U_{p-1}\kron\cdots\kron\U_1\kron\tX\U_0)\vectorize{\G} - \vectorize{\tY}}^2_F\enspace.$$
Let $\M = \U_p\kron \U_{p-1}\kron\cdots\kron\U_1\kron\tX\U_0$. The solution  w.r.t. $\vectorize{\G} $ of this
classical linear least-squares problem is given by $(\M^\top\M)^{-1}\M^\top$. Using the 
mixed-product and inverse properties of the Kronecker product and the column-wise orthogonality of
$\U_1,\cdots,\U_p$ we obtain
$$ \vectorize{\G} = \left(\U_p\kron\cdots\kron\U_1\kron
(\U_0^\top\tX^\top\tX\U_0)^{-1}\U_0^\top\tX^\top\right)\vectorize{\tY}\enspace.\qedhere$$
%The result then follows from observing that $\tX^\top\tX = \X^\top\X+\gamma\I$
%and $\tilde{\ten{Y}}\ttm{1}\tX^\top=\ten{Y}\ttm{1}\X^\top$.
\end{proof}


It follows from Theorem~\ref{thm:core-tensor} that problem~(\ref{eq:holrr-pbm}) 
can be written as 
\begin{align}\label{eq:holrr-pbm3}
\min_{\U_i \in \R^{d_i\times R_i},\atop 0\leq i \leq p} &
\norm{\tY\ttm{1}\mats{\Pi}_0\ttm{2}\cdots\ttm{p+1}\mats{\Pi}_p - \tY}^2_F  \\
\text{\ \ \ \ s.t. }&\U_i^\top\U_i=\I\text{ for } 0\leq i \leq p\nonumber,\\
&\mats{\Pi}_0 = \tX\U_0\left(\U_0\tX^\top\tX\U_0^\top\right)^{-1}\U_0^\top\tX^T,\nonumber\\
&\mats{\Pi}_i = \U_i\U_i^\top \text{ for }1\leq i \leq p\enspace,\nonumber
\end{align}
where we used $ \W\ttm{1}\tX= \tY\ttm{1}\mats{\Pi}_0\ttm{2}\cdots\ttm{p+1}\mats{\Pi}_p$.
Note that $\mats{\Pi}_0$ is the orthogonal projection onto the space spanned by the
columns of $\tX\U_0$ and $\mats{\Pi}_i$ is the orthogonal projection onto the column
space of $\U_i$ for $i\geq 1$.
Hence solving problem~(\ref{eq:holrr-pbm}) is equivalent to finding $p+1$ low-dimensional 
subspaces $U_0,\cdots,U_p$ such that projecting  $\tY$ onto
the spaces $\tX U_0, U_1,\cdots, U_p$ along the corresponding modes is  close to $\tY$.

\paragraph{HOLRR algorithm.}
Since solving problem~(\ref{eq:holrr-pbm3}) for the $p+1$ projections simultaneously is a difficult
non-convex optimization problem we propose to solve it independently for each projection. This 
approach has the benefits of both being computationally efficient and providing good theoretical
approximation guarantees (see Theorem~\ref{thm:approx}). The following proposition gives the
analytic solutions of~(\ref{eq:holrr-pbm3}) when  each projection is considered independently.


\begin{proposition}\label{prop:projection-sol}
For $0\leq i \leq p$, using the definition of $\mats{\Pi}_i$ in~(\ref{eq:holrr-pbm3}), the optimal solution of 
$$\min_{\U_i \in \R^{d_i\times R_i}  }\norm{\tY\ttm{i+1}\mats{\Pi}_i - \tY}^2_F \text{ s.t. } \U_i^\top\U_i = \I$$
is given by the eigenvectors of 
$$
\begin{cases}
(\tX^\top\tX) ^{-1}\tX^\top\tenmat{\tilde{Y}}{1}\tenmat{\tilde{Y}}{1}^\top\tX &\text{ if } i=0\\
\tenmat{\tilde{Y}}{i}\tenmat{\tilde{Y}}{i}^\top &\text{ otherwise}
\end{cases}
$$
that corresponds to the $R_i$ largest eigenvalues.
\end{proposition}
\begin{proof}
For any $0\leq i \leq p$, since $\mats{\Pi}_i$ is a projection we have
$\langle \tY\ttm{1} \mats{\Pi}_i, \tY\rangle = \langle \mats{\Pi}_i\tenmat{\tilde{Y}}{i}, \tenmat{\tilde{Y}}{i}\rangle
= \norm{\mats{\Pi}_i\tenmat{\tilde{Y}}{i}}^2_F$, thus minimizing $\norm{ \tilde{\Y}\ttm{i}\mats{\Pi}_i - \tilde{\Y}}^2_F$
is equivalent to minimizing 
$  \norm{\mats{\Pi}_i\tenmat{\tilde{Y}}{i}}^2_F - 2\langle \mats{\Pi}_i\tenmat{\tilde{Y}}{i}, \tenmat{\tilde{Y}}{i}\rangle
= - \norm{\mats{\Pi}_i\tenmat{\tilde{Y}}{i}}^2_F$.
For $i \geq 1$, we have $\norm{\mats{\Pi}_i\tenmat{\tilde{Y}}{i}}^2_F = \Tr(\U_i^\top \tenmat{\tilde{Y}}{i}\tenmat{\tilde{Y}}{i}^\top \U_i)$
which is maximized by letting the columns of $\U_i$ be the top $R_i$ eigenvectors of
the matrix  $\tenmat{\tilde{Y}}{i}\tenmat{\tilde{Y}}{i}^\top$. 
%
%For $i=0$, let 
%$\tilde{\mats{\Pi}}_0 = \tX\U_0\left(\U_0\tX^\top\tX\U_0^\top\right)^{-1}\U_0^\top\tX^\top$
%and check that $\mats{\Pi}_0\Y\Y^\top\mats{\Pi}_0^\top = 
%\mats{\tilde{\Pi}}_0\tenmat{\tilde{Y}}{1}\tenmat{\tilde{Y}}{1}^\top\tilde{\mats{\Pi}}_0^\top$. 
For $i=0$ we have  
$$
\norm{\mats{\Pi}_0\tenmat{\tilde{Y}}{i}}^2_F = \Tr(\mats{{\Pi}}_0\tenmat{\tilde{Y}}{1}\tenmat{\tilde{Y}}{1}^\top
{\mats{\Pi}}_0^\top)
 =
\Tr\left((\U_0^\top \mat{A} \U_0)^{-1} 
\U_0^\top \mat{B} \U_0\right)
$$
with $\mat{A} = \mat{\tX^\top\tX} $ and 
$\mat{B} =\tX^\top\tenmat{\tilde{Y}}{1}\tenmat{\tilde{Y}}{1}^\top\tX$,
which is maximized by the top $R_0$  eigenvectors of $\mat{A}^{-1}\mat{B}$. 
\end{proof}

The results from Theorem~\ref{thm:core-tensor} and Proposition~\ref{prop:projection-sol} can be
rewritten using the original input matrix $\X$ and output tensor $\Y$.
Since $\tX^\top\tX = \X^\top\X+\gamma\I$ and $\tilde{\ten{Y}}\ttm{1}\tX^\top=\ten{Y}\ttm{1}\X^\top$,
we can rewrite the solution of Theorem~\ref{thm:core-tensor} as 
$$\G = \Y \ttm{1}(\U_0^\top(\X^\top\X+\gamma\I)\U_0)^{-1}\U_0^\top\X^T
\ttm{2}\U_1^\top\ttm{3}\cdots\ttm{p+1}\U_p^\top\enspace.$$ Similarly for
Proposition~\ref{prop:projection-sol} we have
$$(\tX^\top\tX) ^{-1}\tX^\top\tenmat{\tilde{Y}}{1}\tenmat{\tilde{Y}}{1}^\top\tX =
(\X^\top\X + \gamma\I)^{-1}\X^\top\tenmat{Y}{1}\tenmat{Y}{1}^\top\X$$ and one
can check that $\tenmat{\tilde{Y}}{i}\tenmat{\tilde{Y}}{i}^\top = \tenmat{Y}{i}\tenmat{Y}{i}^\top$
for any $i \geq 1$.

The overall Higher-Order Low-Rank Regression procedure (HOLRR)  is summarized in Algorithm~\ref{alg:holrr}.
Note that the Tucker decomposition of the solution returned by HOLRR could be a good initialization point
for an Alternative Least Square method. However, studying the theoretical and experimental properties of
this approach is beyond the scope of this paper and is left for future work.

\begin{algorithm}
   \caption{\texttt{HOLRR}}
   \label{alg:holrr}
\begin{algorithmic}[1]
   \REQUIRE $\mat{X}\in \R^{N \times d_0}$, 
   $\Y \in \R^{N\times d_1\times\cdots\times d_p}$,
   rank $(R_0,R_1,\cdots,R_p)$ and
   regularization parameter $\gamma$.
   \STATE$\U_0\ \leftarrow$ top $R_0$  eigenvectors of\newline
${\hspace{1.5cm} \scriptsize (\X^\top\X + \gamma\I)^{-1}\X^\top\tenmat{Y}{1}\tenmat{Y}{1}^\top\X}$\label{alg:holrr-U0} 
   \FOR{$i=1$ {\bfseries to} $p$}
   \STATE$\U_i \leftarrow$ top $R_i$ eigenvectors of $\tenmat{Y}{i}\tenmat{Y}{i}^\top$\label{alg:holrr-Ui} 
   \ENDFOR
   \STATE $\mat{M} = \left(\U_0^\top(\X^\top\X+\gamma\I)\U_0\right)^{-1}\U_0^\top\X^\top$
	\STATE $\G \leftarrow \Y \ttm{1}\mat{M}\ttm{2}\U_1^\top\ttm{3}\cdots\ttm{p+1}\U_p^\top$
	\STATE  {\bfseries return}  $\G\ttm{1}\U_0\ttm{2}\cdots\ttm{p+1}\U_p$
\end{algorithmic}
\end{algorithm}

 
conclusion.tex
\section{Conclusion}
\label{sec:ccl}
We proposed a low-rank multilinear regression model for tensor-structured output data.
We developed a fast and efficient algorithm to tackle the multilinear rank penalized
minimization problem and provided theoretical approximation guarantees. Experimental
results showed that capturing low-rank structure in the output data can help to 
improve tensor regression performance.
guarantees.tex
%\section{Theoretical Guarantees}


\subsection{Theoretical Analysis}

\paragraph{Complexity analysis.} We compare the computational complexity
of the LRR and  HOLRR algorithms.
% with the ones of RLS and LRR (after vectorization of the outputs). 
%The complexity of RLS is $\bigo{N(d_0)^2 + (d_0)^3}$.
For both algorithms the computational cost for matrix multiplications is asymptotically the same
(dominated by the products $\X^\top\X$ and $\tenmat{Y}{1}\tenmat{Y}{1}^\top$). 
We thus focus our analysis on the computational cost of
matrix inversions and truncated singular value decompositions. For the two methods
the inversion of the matrix $\X^\top\X +\gamma\I$ can be done in $\bigo{(d_0)^3}$. The LRR
method needs to compute the $R$ dominant eigenvectors of a matrix of
the same size as $\tenmat{Y}{1}^\top \tenmat{Y}{1}$ (see Section~\ref{sec:LRR})
which can be done in $\bigo{R(d_1d_2\cdots d_p)^2}$. The HOLRR algorithm needs
to compute the $R_i$ dominant eigenvectors of the matrix $\tenmat{Y}{i} \tenmat{Y}{i}^\top$
for $i\in [p]$ and to invert the extra matrix $\U_0^\top(\X^\top\X+\gamma\I)\U_0$ of
size $R_0\times R_0$ leading to a complexity of $\bigo{(R_0)^3 + \max_i\{R_id_i^2\}}$. 
We can conclude that 
%while the HOLRR algorithm is only lightly more computationally expensive than RLS
HOLRR is far more efficient than LRR when the output dimensions $d_1,\cdots,d_p$ are large.

\paragraph{Approximation guarantees.}
The following theorem shows that the HOLRR algorithm proposed in the previous section
is a $(p+1)$-approximation algorithm for problem~(\ref{eq:holrr-pbm}).
\begin{theorem}
\label{thm:approx}
Let $\W^*$ be a solution of problem~(\ref{eq:holrr-pbm}) and
let $\W$ be the regression tensor  returned by Algorithm~\ref{alg:holrr}.
If $\mathcal{L}:\R^{d_0\times \cdots \times d_p } \to \R$ denotes the objective 
function of~(\ref{eq:holrr-pbm}) w.r.t. $\W$ then
$\mathcal{L}(\W) \leq  (p+1) \mathcal{L}(\W^*) .$
\end{theorem}
\begin{proof}
Let $\U_0,\cdots,\U_p$ be the matrices defined in Algorithm~\ref{alg:holrr} and
let $\mats{\Pi}_0,\cdots,\mats{\Pi}_p$ be the projection matrices defined in
problem~(\ref{eq:holrr-pbm3}).
We have

\begin{align*}
\mathcal{L}(\W) &= \norm{\W\ttm{1}\tilde{\X} - \tilde{\Y}}^2_F\\
&= \norm{\tY\ttm{1}\mats{\Pi}_0\ttm{2}\cdots\ttm{p+1}\mats{\Pi}_p - \tY}^2_F \\
&\leq \norm{\tY\ttm{1}\mats{\Pi}_0\ttm{2}\cdots\ttm{p+1}\mats{\Pi}_p 
                     - \tY\ttm{1}\mats{\Pi}_0\ttm{2}\cdots\ttm{p}\mats{\Pi}_{p-1}}^2_F \\
&+ \norm{\tY\ttm{1}\mats{\Pi}_0\ttm{2}\cdots\ttm{p}\mats{\Pi}_{p-1}
               - \tY\ttm{1}\mats{\Pi}_0\ttm{2}\cdots\ttm{p-1}\mats{\Pi}_{p-2}}^2_F \\
 &+\cdots 
 + \norm{\tY\ttm{1}\mats{\Pi}_0 - \tY}^2_F\\
&\leq \norm{\tY\ttm{p+1}\mats{\Pi}_p  - \tY}^2_F + \norm{\tY\ttm{p}\mats{\Pi}_{p-1}-\tY}^2_F
+ \cdots +   \norm{\tY\ttm{1}\mats{\Pi}_0 - \tY}^2_F\enspace,
\end{align*}
%
where we used the fact that $\norm{\mats{\Pi}\M}_F \leq \norm{\M}_F$ for any projection matrix $\mats{\Pi}$.
Furthermore, it follows from Theorem~\ref{thm:core-tensor} that 
$\W^* = \Y \ttm{1}\mats{\Pi}_0^* \ttm{2} \cdots \ttm{p+1}\mats{\Pi}_p^*$
for some projection matrices $\mats{\Pi}_i^*$ satisfying the constraints from problem~(\ref{eq:holrr-pbm3}).
By Proposition~\ref{prop:projection-sol} each summand in the last inequality is minimal w.r.t. $\mats{\Pi}_i$, hence
\begin{align*}
\norm{\tY\ttm{i}\mats{\Pi}_i - \tY}^2_F &\leq \norm{\tY\ttm{i}\mats{\Pi}^*_i - \tY}^2_F\\
&\leq  \norm{\tY\ttm{1}\mats{\Pi}^*_0\ttm{2}\cdots \ttm{p+1}\mats{\Pi}^*_p - \tY}^2_F \\
 &= \norm{\W^*\ttm{1}\tX- \tY}^2_F \\
 &= \mathcal{L}(\W^*)
\end{align*}
which concludes the proof.
\end{proof}

One direct consequence of this theorem is that if there exists a $\W^*$ such
that $\W^*\ttm{1}\X = \Y$ then Algorithm~\ref{alg:holrr} (using $\gamma=0$)
will return a tensor $\W$ satisfying $\W\ttm{1}\X = \Y$.

%\subsection{Rademacher Complexity}
%\input{rademacher}




intro.tex
\section{Introduction}


Recently, there has been an increasing interest in adapting machine learning and statistical methods to tensors.
 Data with a natural tensor structure is encountered in many scientific areas
 including neuroimaging~\citep{Zhou2013}, signal processing~\citep{Cichocki2009}, spatio-temporal analysis~\citep{Bahadori2014} 
 and computer vision~\citep{Lu2013}.
 %
Extending multivariate regression methods to tensors is one of the challenging task in this area.
Most existing works extend linear models to the multilinear setting and focus on the tensor structure 
of the input data (e.g.~\citet{signoretto2013a}). Little has been done however to investigate  learning methods 
for \emph{tensor-structured output data}. 

We consider a multilinear regression task where inputs are vectors and outputs are tensors. 
In order to leverage the tensor structure of the output data, we formulate the problem as the
minimization of a least squares criterion 
subject to a \emph{multilinear rank} constraint on the regression tensor. The rank constraint enforces 
the model to capture the low-rank structure of the outputs and to explain
the dependencies between inputs and outputs in a low-dimensional multilinear subspace. 

Unlike previous work we do not use a convex relaxation of this difficult non-convex optimization problem. Instead we  design
an efficient approximation algorithm (HOLRR) for which we are able to provide good approximation guarantees.
We also present a kernelized version of HOLRR which extends our model to the nonlinear
setting. Experiments on synthetic and real data shows that HOLRR obtains better predictive accuracy
while being computationally very competitive. We also present an image recovery experiment
which gives an illustrative insight on the effects of multilinear rank regularization.

%experiments showing that HOLRR   is computationally very competitive.


\paragraph{Related work.}  In the context of multi-task learning, \citet{romera2013} have proposed a linear model using a 
tensor-rank penalization of a least squares criterion  to take into account the multi-modal interactions between tasks. Their
approach relies on a convex relaxation of the multlinear rank constraint using the trace norms of the matricizations. They also
propose a non-convex approach (MLMT-NC) but it is computationally very expensive.
%
\citet{Bahadori2014} have proposed a greedy algorithm to solve a low-rank tensor learning problem in the context of multivariate 
spatio-temporal data analysis. The linear model they assume is different from the one we propose, it is specifically 
designed for spatio-temporal data and does not fit into the general tensor-valued regression framework we consider here.



\paragraph{Paper outline.}
We provide some background on low-rank regression and tensors in Section~\ref{sec:prelim}. 
In Section~\ref{sec:algo} we introduce the tensor response regression problem,
we develop the HOLRR algorithm to tackle the minimization problem with mulitlinear rank constraint, and
we prove that this algorithm has good approximation guarantees. A kernelized version of HOLRR is
provided in Section~\ref{sec:kernel}. Finally, we assess the performances of our method through simulation
study and real data analysis in Section~\ref{sec:xp}.

\comment{
Multiple output regression is the task of predicting $p$ multiple response variables from a shared set of $d$ predictor variables. 
It is well known that the ordinary least squares method in this context is equivalent to performing $p$ independent linear regressions. 
To address this issue the reduced-rank regression model learns a low-rank approximation of the regression matrix that maps 
the covariate vectors in $\mathbb{R}^d$ to the response vectors in $\mathbb{R}^p$ 
\citep{Izenman75,Reinsel1998,Vounou2010,Chen2012,Zheng2014}. 
Rank regularization has two main advantages: it reduces the number of model parameters, and it takes advantage of the linear 
dependencies between the response variables. 
To solve a regression task with tensor-structured output, one could vectorize the outputs and perform a standard reduced-rank 
regression. However we argue that this approach would fail to capture the higher level dependencies expected to occur 
in tensor-structured outputs.

In this paper, we show how the reduced-rank model can be extended to tensor responses by minimizing a least squares criterion 
subject to a multilinear rank constraint on the regression tensor. 
The rank constraint enforces the model to capture the low-rank structure of the outputs and to explain
the dependencies between inputs and outputs in a low-dimensional multilinear space. 

\paragraph{Related work.}  In the context of multi-task learning, \citet{romera2013} have proposed a linear model using a 
tensor-rank penalization of a least squares criterion  to take into account the multi-modal interactions between tasks. 
\citet{Bahadori2014} have proposed a greedy algorithm to solve a low-rank tensor learning problem in the context of multivariate 
spatio-temporal data analysis. The linear model they assume is different from the one we propose, it is specifically 
designed for spatio-temporal data and does not fit into the general tensor-valued regression framework we consider here.
Both approaches rely on a convex relaxation of the multlinear rank constraint using the trace norms of the matricizations. 
 In contrast, we directly tackle the tensor rank constraint and we are able to provide good theoretical approximation
 guarantees. \citet{Foygel2012} have introduced a nonparametric extension of the linear 
 reduced-rank framework, and a kernel extension of reduced-rank regression was studied in~\citet{Mukherjee2011}. 
 However, only the case of vector outputs has been considered. 


\paragraph{Paper outline.}
We provide some background on reduced-rank regression and tensors in Section~\ref{sec:prelim}. 
In Section~\ref{sec:algo} we present the extension of reduced-rank regression for tensor responses,
we develop the HOLRR algorithm to tackle the minimization problem with mulitlinear rank constraint and
we prove that this algorithm has good approximation guarantees. In Section~\ref{sec:kernel}, we provide
a kernelized version of HOLRR. We assess the performances of our method through simulation
study and real data analysis in Section~\ref{sec:xp}. Finally, Section~\ref{sec:cccl}  concludes 
the paper and suggests directions for future work.

}
\comment{

On the other hand, almost all the efforts in reduced-rank regression have focused on linear models. 
Non-parametric reduced-rank regression has received little interest, even in the case of vector responses. 
In this paper, we propose to handle tensor-structured responses by relying on the operator-valued kernel framework. Operator-valued kernels extend scalar-valued kernels to multi-output situations. Since the seminal work of~\citet{Micchelli-2005a}, operator-valued kernels have been successfully used for various problems in statistics and machine learning, such as multitask learning~\citep{Evgeniou-2005}, multivariate regression~\citep{Sindhwani13}, functional regression~\citep{Kadri-2010}, structured output prediction~\citep{Brouard11,Kadri13} and network inference~\citep{Lim2014}. However, to the best of our knowledge, using operator-valued kernels to deal with higher order tensor-structured outputs has not been investigated yet. We show how operator-valued kernels can be a powerful tool for the design of nonparametric tensor regression procedures by introducing the notion of tensor-valued kernels. 

\medskip

\textbf{Related Work.}
\medskip

\textbf{Paper Outline.} %The paper is organized as follows. 
We provide some background on reduced-rank regression, kernels and tensors in Section~\ref{sec:BN}. In Section~\ref{sec:TRR}, we show how the reduced-rank regression model can be extended to tensor-structured output and we introduce the notion of tensor-valued kernels  which we use to extend our model to the nonparametric setting. In Section~\ref{sec:learning}, we propose two algorithms to learn the tensor-valued regression function from data. The first one attempts to solve the tensor reduced-rank problem for identity-based tensor-valued kernels, while the second jointly learns the tensor regressor and the tensor-valued kernel. 
%Finally, we provide an empirical evaluation of the performance of our learning algorithms which demonstrates their effectiveness on simulated and real data in Section~\ref{sec:xp}, followed by a brief conclusion.
We provide an empirical evaluation of the performance of our learning algorithms which demonstrates their effectiveness on simulated and real data  in Section~\ref{sec:xp}. Finally, Section~\ref{sec:conclusion}  concludes the paper and suggests directions for future work.
}


kernel.tex
\section{HOLRR Kernel Extension}
\label{sec:kernel}
In this section we provide a kernelized version of the HOLRR algorithm. We proceed by analyzing
how the algorithm would be instantiated in a feature space and we show that all the 
steps involved can be performed using the Gram matrix of the input data without
having to explicitly compute the feature map. 

Let $\phi: \R^{d_0} \to \R^L$ be a feature map and let $\mats{\Phi} \in \R^{n\times L}$ be the
matrix with lines $\phi(\vec{x}\nsample{n})^\top$ for $n\in[N]$. The higher-order low-rank 
regression problem in the feature space boils down to the minimization problem
\begin{align}\label{eq:kholrr-pbm}
\MoveEqLeft\min_{\W\in \R^{L\times d_1 \times \cdots \times d_p}} \norm{\W\ttm{1}\mats{\Phi} - \Y}^2_F 
+ \gamma \norm{\W}^2_F \\
&\text{s.t. }\rank(\W) \leq  (R_0, R_1, \cdots, R_p)\enspace.\nonumber
\end{align}

Following the HOLRR algorithm proposed in Section~\ref{sec:holrr} one needs to compute
the top $R_0$ eigenvectors of the $L\times L$ matrix
$( \mats{\Phi}^\top\mats{\Phi}  + \gamma\I)^{-1}\mats{\Phi}^\top\tenmat{Y}{1}\tenmat{Y}{1}^\top\mats{\Phi} $.
The following proposition shows that this can be done using the Gram matrix $\K = \mats{\Phi} \mats{\Phi} ^\top$
without having to explicitly  compute the feature map $\phi$.

\begin{proposition}
If $\vecs{\alpha}\in\R^N$ is an eigenvector with eigenvalue $\lambda$ of the matrix 
\begin{equation}
\label{eq:gram-eigvec}
(\K+\gamma\I)^{-1}\tenmat{Y}{1}\tenmat{Y}{1}^\top \K\enspace,
\end{equation}
then $\vec{v} = \mats{\Phi}^\top \vecs{\alpha} \in\R^L$ is an eigenvector with eigenvalue $\lambda$ of
the matrix
$( \mats{\Phi}^\top\mats{\Phi}  + \gamma\I)^{-1}\mats{\Phi}^\top\tenmat{Y}{1}\tenmat{Y}{1}^\top\mats{\Phi} .$
\end{proposition}
\begin{proof}
Let  $\vecs{\alpha}\in\R^N$ be the eigenvector from the hypothesis. We have
\begin{align*}
\lambda\vec{v}
&=
\mats{\Phi}^\top (\lambda \vecs{\alpha}) 
=
\mats{\Phi}^\top\left( (\K+\gamma\I)^{-1}\tenmat{Y}{1}\tenmat{Y}{1}^\top \K \right)\vecs{\alpha}\\
&=
\mats{\Phi}^\top( \mats{\Phi} \mats{\Phi}^\top + \gamma\I)^{-1}
\tenmat{Y}{1}\tenmat{Y}{1}^\top\mats{\Phi} \mats{\Phi}^\top \vecs{\alpha}\\
&=
\left(( \mats{\Phi}^\top\mats{\Phi}  + \gamma\I)^{-1}
\mats{\Phi}^\top\tenmat{Y}{1}\tenmat{Y}{1}^\top\mats{\Phi} \right)\vec{v} \enspace.\qedhere
\end{align*}
\end{proof}

Let $\mat{A}$ be the top $R_0$ eigenvectors of the matrix~(\ref{eq:gram-eigvec}).
When working with the feature map $\phi$, it follows from the previous proposition that
 line~1 in Algorithm~\ref{alg:holrr} is  equivalent
to choosing $\U_0 = \mats{\Phi}^\top \mat{A}\in \R^{L\times R_0}$, while the updates
in line~3 stay the same. The regression tensor $\W\in\R^{L\times d_1\times\cdots\times d_p}$
 returned by this algorithm  is then equal to
\begin{align*}
\W = \Y\ttm{1}  \mat{P} \ttm{2} \U_1\U_1^\top \ttm{2} \cdots \ttm{p+1} \U_p\U_p^\top\enspace,
\end{align*}
where 
\begin{align*}
\mat{P} 
&= 
\mats{\Phi}^\top\mat{A}\left(\mat{A}^\top\mats{\Phi}
(\mats{\Phi}^\top\mats{\Phi}+\gamma\I) \mats{\Phi}^\top\mat{A}\right)^{-1}\mat{A}^\top\mats{\Phi}\mats{\Phi}^\top\\
&=
\mats{\Phi}^\top\mat{A}\left(\mat{A}^\top\mats{\Phi}\mats{\Phi}^\top
(\mats{\Phi} \mats{\Phi}^\top+\gamma\I)\mat{A}\right)^{-1}\mat{A}^\top\mats{\Phi}\mats{\Phi}^\top\\
&=
\mats{\Phi}^\top\mat{A}\left(\mat{A}^\top\K
(\K+\gamma\I)\mat{A}\right)^{-1}\mat{A}^\top\K\enspace.
\end{align*}

Suppose now that the feature map $\phi$ is induced by a kernel $k:\R^{d_0}\times\R^{d_0}\to \R$.
The prediction for an input vector $\vec{x}$ is then given by 
$\W\ttv{1}\vec{x} = \ten{C}\ttv{1}\vec{k}_{\vec{x}}$ where
the $n$th component of $\vec{k}_{\vec{x}}\in \R^N$ is $k(\vec{x}_n,\vec{x})$ and
the tensor $\ten{C}\in\R^{N\times  d_1\times \cdots \times d_p}$ is defined by
%
$$\ten{C} = \G \ttm{1}\mat{A}\ttm{2}\U_1\ttm{2}\cdots\ttm{p+1}\U_p^\top\enspace,$$ 
%
with
$\G = \Y\ttm{1} \left(\mat{A}^\top\K(\K+\gamma\I)\mat{A}\right)^{-1}\mat{A}^\top\K
\ttm{2}\U_2^\top\ttm{3}\cdots\ttm{p+1}\U_p$. 
%Note that the size of the core tensor
%$\G$ is $R_0\times R_1\times \cdots \times R_p$ as in the linear case but the first
%factor  $\mat{A}$ is now of size $N\times R_0$.

Now, let $\mathcal{H}$ be the reproducing kernel Hilbert space associated with the kernel $k$.
The overall procedure for kernelized HOLRR is summarized in Algorithm~\ref{alg:kholrr}. This
algorithm returns the tensor $\ten{C} \in \R^{N\times d_1 \times \cdots \times d_p}$ defining 
the regression function 
%
%\vspace{-0.5cm}%
%
$$f: \vec{x} \mapsto \ten{C}\ttv{1}\vec{k}_{\vec{x}} = \sum_{n=1}^N k(\vec{x},\vec{x}^{(n)}) \ten{C}^{(n)}\enspace,$$
%
%\vspace{-0.5cm}
%
where $\ten{C}^{(n)} = \ten{C}_{n:\cdots:}\in \R^{d_1\times\cdots\times d_p}$.

\begin{algorithm}
   \caption{\texttt{Kernelized HOLRR}}
   \label{alg:kholrr}
\begin{algorithmic}[1]
   \REQUIRE Gram matrix $\mat{K}\in \R^{N \times N}$, 
   $\Y \in \R^{N\times d_1\times\cdots\times d_p}$,
   rank $(R_0,R_1,\cdots,R_p)$ and
   regularization parameter $\gamma$.
   \STATE\label{alg:holrr-U0} $\mat{A}\ \leftarrow$ top $R_0$  eigenvectors of\newline
${\hspace{1.5cm} \scriptsize (\K+\gamma\I)^{-1}\tenmat{Y}{1}\tenmat{Y}{1}^\top \K}$
   \FOR{$i=1$ {\bfseries to} $p$}
   \STATE\label{alg:holrr-U0} $\U_i \leftarrow$ top $R_i$ eigenvectors of $\tenmat{Y}{i}\tenmat{Y}{i}^\top$
   \ENDFOR
   \STATE $\mat{M}\leftarrow \left(\mat{A}^\top\K(\K+\gamma\I)\mat{A}\right)^{-1}\mat{A}^\top\K$
	\STATE $\G \leftarrow \Y \ttm{1}\mat{M} \ttm{2}\U_1^\top\ttm{3}\cdots\ttm{p+1}\U_p^\top$
	\STATE  {\bfseries return}  $\ten{C} = \G\ttm{1}\mat{A}\ttm{2}\U_1\ttm{3}\cdots\ttm{p+1}\U_p$
\end{algorithmic}
\end{algorithm}
macros.tex

%% Macros to make visible comments while writing the paper
\newcommand{\hachem}[1]{\todo[inline,color=green!40]{{\it Hachem:~}#1}}
\newcommand{\guillaume}[1]{\todo[inline,color=blue!40]{{\it Guillaume:~}#1}}
\renewcommand{\comment}[1]{}

%% Mathematical sets with names
\newcommand{\R}{\mathbb{R}}

%% Probability notations
\newcommand{\Prob}{\mathbb{P}}
\newcommand{\Esp}{\mathbb{E}}
\newcommand{\Var}{\mathbb{V}}

%% Vectors, matrices and tensors
\renewcommand{\vec}[1]{\ensuremath{\mathbf{#1}}}
\newcommand{\vecs}[1]{\ensuremath{\mathbf{\boldsymbol{#1}}}}
\newcommand{\mat}[1]{\ensuremath{\mathbf{#1}}}
\newcommand{\mats}[1]{\ensuremath{\mathbf{\boldsymbol{#1}}}}
\newcommand{\ten}[1]{\mat{\ensuremath{\boldsymbol{\mathcal{#1}}}}}
% tensor times matrix
\newcommand{\ttm}[1]{\times_{#1}}
% tensor times vector
\newcommand{\ttv}[1]{\bullet_{#1}}
% matricisation
\newcommand{\tenmat}[2]{\mat{#1}_{(#2)}}
\newcommand{\tenmatpar}[2]{(#1)_{(#2)}}

% sample number
\newcommand{\nsample}[1]{^{(#1)}}

% common matrices and tensors
\newcommand{\T}{\ten{T}}
\newcommand{\W}{\ten{W}}
\newcommand{\G}{\ten{G}}
\newcommand{\Y}{\ten{Y}}
\newcommand{\U}{\mat{U}}
\newcommand{\X}{\mat{X}}
\newcommand{\I}{\mat{I}}
\newcommand{\M}{\mat{M}}
\newcommand{\K}{\mat{K}}
\renewcommand{\v}{\vec{v}}
\newcommand{\x}{\vec{x}}
% adding a tilde for input matrix and output tensor
\newcommand{\tX}{\tilde{\X}}
\newcommand{\tY}{\tilde{\Y}}



%% Math environments
\newtheorem*{theorem*}{Theorem}
\newtheorem{theorem}{Theorem}
\newtheorem{proposition}{Proposition}
\newtheorem{lemma}{Lemma}
\newtheorem{definition}{Definition}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{condition}{Condition}
\newtheorem{claim}{Claim}

%% Miscellaneous math stuff
\DeclareMathOperator*{\Tr}{Tr} %trace
\DeclareMathOperator*{\argmax}{arg\,max}
\newcommand{\vectorize}[1]{\mathrm{vec}(#1)}
\DeclareMathOperator*{\reshape}{reshape}
\newcommand{\eqdef}{\triangleq}
\newcommand{\vtimes}[1]{\ \overline{\times}_{#1}\ }
\DeclareMathOperator*{\rank}{rank}
\newcommand{\norm}[1]{\|#1\|}
\newcommand{\kron}{\otimes}
\newcommand{\Ker}{\mathrm{Ker}}
\renewcommand{\Im}{\mathrm{Im}}
\newcommand{\bigo}[1]{\mathcal{O}(#1)}
\newcommand{\singv}{\mathfrak{s}}
main.tex

% fix ordering of figure* bug
%\RequirePackage{fixltx2e}

\documentclass{article}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ICML default packages
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% use Times
\usepackage{times}
% For figures
\usepackage{graphicx} % more modern
%\usepackage{epsfig} % less modern
\usepackage{subfigure} 

% For citations
\usepackage{natbib}

% For algorithms
\usepackage{algorithm}
\usepackage{algorithmic}


\usepackage{hyperref}

\newcommand{\theHalgorithm}{\arabic{algorithm}}

\makeatletter
\providecommand\theHALG@line{\thealgorithm.\arabic{ALG@line}}
\makeatother

%\usepackage{icml2016} 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% our packages
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{mathtools}
\usepackage{lmodern}

% fix input in algorithm environment
\renewcommand{\algorithmicrequire}{\textbf{Input:}}

% for notes (example: \hachem{We need to change all the notations...} \guillaume{No way!})
\usepackage[textsize=footnotesize,color=green!40]{todonotes}

% for \llbracket and \rrbracket
\usepackage{stmaryrd}

% avoid figure taking the whole page
\renewcommand{\floatpagefraction}{.95}%
%\usepackage{stfloats}

% for nice looking hyperlinks
\usepackage{xcolor}
\hypersetup{
    colorlinks,
    linkcolor={red!80!black},
    citecolor={blue!80!black},
    urlcolor={blue!80!black}
}


\newcommand{\comment}[1]{}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% our macros
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\input{macros.tex}

% The \icmltitle you define below is probably too long as a header.
% Therefore, a short form for the running title is supplied here:
\title{Higher-Order Low-Rank Regression}
\author{Guillaume Rabusseau\footnote{Contact author: \texttt{guillaume.rabusseau@lif.univ-mrs.fr}.} \qquad Hachem Kadri \\ Aix-Marseille University  }


\begin{document}
\maketitle


\begin{abstract} 
\input{abstract.tex}
\end{abstract} 

\input{intro.tex}

%\input{notation.tex}
\input{prelim.tex}
\input{algo.tex}
\input{guarantees.tex}
\input{kernel.tex}
\input{xp.tex}
\input{conclusion.tex}

\newpage

\newpage
\bibliographystyle{icml2016}
\bibliography{bib_TensorKernels}

%\section{Supplementary Material}
%\twocolumn[
%\input{supmat.tex}
%]
\end{document}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ICML guidelines
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document} 

\twocolumn[
\icmltitle{Submission and Formatting Instructions for \\ 
           International Conference on Machine Learning (ICML 2016)}

% It is OKAY to include author information, even for blind
% submissions: the style file will automatically remove it for you
% unless you've provided the [accepted] option to the icml2016
% package.
\icmlauthor{Your Name}{email@yourdomain.edu}
\icmladdress{Your Fantastic Institute,
            314159 Pi St., Palo Alto, CA 94306 USA}
\icmlauthor{Your CoAuthor's Name}{email@coauthordomain.edu}
\icmladdress{Their Fantastic Institute,
            27182 Exp St., Toronto, ON M6H 2T1 CANADA}

% You may provide any keywords that you 
% find helpful for describing your paper; these are used to populate 
% the "keywords" metadata in the PDF but will not be shown in the document
\icmlkeywords{boring formatting information, machine learning, ICML}

\vskip 0.3in
]

\begin{abstract} 
The purpose of this document is to provide both the basic paper template and
submission guidelines. Abstracts should be a single paragraph, between 4--6 sentences long, ideally.  Gross violations will trigger corrections at the camera-ready phase.
\end{abstract} 

\section{Electronic Submission}
\label{submission}

Submission to ICML 2016 will be entirely electronic, via a web site
(not email).  Information about the submission process and \LaTeX\ templates
are available on the conference web site at:
\begin{center}
\textbf{\texttt{http://icml.cc/2016/}}
\end{center}
Send questions about submission and electronic templates to
\texttt{icml2016pc@gmail.com}.

The guidelines below will be enforced for initial submissions and
camera-ready copies.  Here is a brief summary:
\begin{itemize}
\item Submissions must be in PDF.
\item The maximum paper length is \textbf{8 pages excluding references and acknowledgements, and 10 pages
  including references and acknowledgements} (pages 9 and 10 must contain only references and acknowledgements).
\item Do \textbf{not include author information or acknowledgements} in your initial
submission. 
\item Your paper should be in \textbf{10 point Times font}.
\item Make sure your PDF file only uses Type-1 fonts.
\item Place figure captions {\em under} the figure (and omit titles from inside
the graphic file itself).  Place table captions {\em over} the table.
\item References must include page numbers whenever possible and be as complete
as possible.  Place multiple citations in chronological order.  
\item Do not alter the style template; in particular, do not compress the paper
format by reducing the vertical spaces.
\item Keep your abstract brief and self-contained, one
   paragraph and roughly 4--6 sentences.  Gross violations will require correction at the camera-ready phase.
  Title should have content words capitalized.
  

\end{itemize}

\subsection{Submitting Papers}

{\bf Paper Deadline:} The deadline for paper submission to ICML 2016
is at \textbf{23:59 Universal Time (3:59 p.m.\ Pacific Standard Time) on February 5, 2016}.
If your full submission does not reach us by this time, it will 
not be considered for publication. There is no separate abstract submission.

{\bf Anonymous Submission:} To facilitate blind review, no identifying
author information should appear on the title page or in the paper
itself.  Section~\ref{author info} will explain the details of how to
format this.

{\bf Simultaneous Submission:} ICML will not accept any paper which,
at the time of submission, is under review for another conference or
has already been published. This policy also applies to papers that
overlap substantially in technical content with conference papers
under review or previously published. ICML submissions must not be
submitted to other conferences during ICML's review period. Authors
may submit to ICML substantially different versions of journal papers
that are currently under review by the journal, but not yet accepted
at the time of submission. Informal publications, such as technical
reports or papers in workshop proceedings which do not appear in
print, do not fall under these restrictions.

\medskip

To ensure our ability to print submissions, authors must provide their
manuscripts in \textbf{PDF} format.  Furthermore, please make sure
that files contain only Type-1 fonts (e.g.,~using the program {\tt
  pdffonts} in linux or using File/DocumentProperties/Fonts in
Acrobat).  Other fonts (like Type-3) might come from graphics files
imported into the document.

Authors using \textbf{Word} must convert their document to PDF.  Most
of the latest versions of Word have the facility to do this
automatically.  Submissions will not be accepted in Word format or any
format other than PDF. Really. We're not joking. Don't send Word.

Those who use \textbf{\LaTeX} to format their accepted papers need to pay close
attention to the typefaces used.  Specifically, when producing the PDF by first
converting the dvi output of \LaTeX\ to Postscript the default behavior is to
use non-scalable Type-3 PostScript bitmap fonts to represent the standard
\LaTeX\ fonts. The resulting document is difficult to read in electronic form;
the type appears fuzzy. To avoid this problem, dvips must be instructed to use
an alternative font map.  This can be achieved with the following two commands:

{\footnotesize
\begin{verbatim}
dvips -Ppdf -tletter -G0 -o paper.ps paper.dvi
ps2pdf paper.ps
\end{verbatim}}
Note that it is a zero following the ``-G''.  This tells dvips to use
the config.pdf file (and this file refers to a better font mapping).

A better alternative is to use the \textbf{pdflatex} program instead of
straight \LaTeX. This program avoids the Type-3 font problem, however you must
ensure that all of the fonts are embedded (use {\tt pdffonts}). If they are
not, you need to configure pdflatex to use a font map file that specifies that
the fonts be embedded. Also you should ensure that images are not downsampled
or otherwise compressed in a lossy way.

Note that the 2016 style files use the {\tt hyperref} package to
make clickable links in documents.  If this causes problems for you,
add {\tt nohyperref} as one of the options to the {\tt icml2016}
usepackage statement.

\subsection{Reacting to Reviews}

We will continue the ICML tradition in which the authors are given the
option of providing a short reaction to the initial reviews. These
reactions will be taken into account in the discussion among the
reviewers and area chairs.

\subsection{Submitting Final Camera-Ready Copy}

The final versions of papers accepted for publication should follow the
same format and naming convention as initial submissions, except of
course that the normal author information (names and affiliations)
should be given.  See Section~\ref{final author} for details of how to
format this.

The footnote, ``Preliminary work.  Under review by the International
Conference on Machine Learning (ICML).  Do not distribute.'' must be
modified to ``\textit{Proceedings of the
$\mathit{33}^{rd}$ International Conference on Machine Learning},
New York, NY, USA, 2016.  JMLR: W\&CP volume 48. 
Copyright 2016 by the author(s).''

For those using the \textbf{\LaTeX} style file, simply change
$\mathtt{\backslash usepackage\{icml2016\}}$ to 

$$\mathtt{\backslash usepackage[accepted]\{icml2016\}}$$

\noindent
Authors using \textbf{Word} must edit the
footnote on the first page of the document themselves.

Camera-ready copies should have the title of the paper as running head
on each page except the first one.  The running title consists of a
single line centered above a horizontal rule which is $1$ point thick.
The running head should be centered, bold and in $9$ point type.  The
rule should be $10$ points above the main text.  For those using the
\textbf{\LaTeX} style file, the original title is automatically set as running
head using the {\tt fancyhdr} package which is included in the ICML
2016 style file package.  In case that the original title exceeds the
size restrictions, a shorter form can be supplied by using

\verb|\icmltitlerunning{...}|

just before $\mathtt{\backslash begin\{document\}}$.
Authors using \textbf{Word} must edit the header of the document themselves.

\section{Format of the Paper} 
 
All submissions must follow the same format to ensure the printer can
reproduce them without problems and to let readers more easily find
the information that they desire.

\subsection{Length and Dimensions}

Papers must not exceed eight (8) pages, including all figures, tables,
and appendices, but excluding references and acknowledgements. When references and acknowledgements are included,
the paper must not exceed ten (10) pages.
Acknowledgements should be limited to grants and people who contributed to the paper.
Any submission that exceeds 
this page limit or that diverges significantly from the format specified 
herein will be rejected without review.

The text of the paper should be formatted in two columns, with an
overall width of 6.75 inches, height of 9.0 inches, and 0.25 inches
between the columns. The left margin should be 0.75 inches and the top
margin 1.0 inch (2.54~cm). The right and bottom margins will depend on
whether you print on US letter or A4 paper, but all final versions
must be produced for US letter size.

The paper body should be set in 10~point type with a vertical spacing
of 11~points. Please use Times typeface throughout the text.

\subsection{Title}

The paper title should be set in 14~point bold type and centered
between two horizontal rules that are 1~point thick, with 1.0~inch
between the top rule and the top edge of the page. Capitalize the
first letter of content words and put the rest of the title in lower
case.

\subsection{Author Information for Submission}
\label{author info}

To facilitate blind review, author information must not appear.  If
you are using \LaTeX\/ and the \texttt{icml2016.sty} file, you may use
\verb+\icmlauthor{...}+ to specify authors.  The author information
will simply not be printed until {\tt accepted} is an argument to the
style file. Submissions that include the author information will not
be reviewed.

\subsubsection{Self-Citations}

If your are citing published papers for which you are an author, refer
to yourself in the third person. In particular, do not use phrases
that reveal your identity (e.g., ``in previous work \cite{langley00}, we 
have shown \ldots'').

Do not anonymize citations in the reference section by removing or
blacking out author names. The only exception are manuscripts that are
not yet published (e.g. under submission). If you choose to refer to
such unpublished manuscripts \cite{anonymous}, anonymized copies have 
to be submitted
as Supplementary Material via CMT. However, keep in mind that an ICML
paper should be self contained and should contain sufficient detail
for the reviewers to evaluate the work. In particular, reviewers are
not required to look a the Supplementary Material when writing their
review.

\subsubsection{Camera-Ready Author Information}
\label{final author}

If a paper is accepted, a final camera-ready copy must be prepared.
%
For camera-ready papers, author information should start 0.3~inches
below the bottom rule surrounding the title. The authors' names should
appear in 10~point bold type, electronic mail addresses in 10~point
small capitals, and physical addresses in ordinary 10~point type.
Each author's name should be flush left, whereas the email address
should be flush right on the same line. The author's physical address
should appear flush left on the ensuing line, on a single line if
possible. If successive authors have the same affiliation, then give
their physical address only once.

A sample file (in PDF) with author names is included in the ICML2016 
style file package.

\subsection{Abstract}

The paper abstract should begin in the left column, 0.4~inches below
the final address. The heading `Abstract' should be centered, bold,
and in 11~point type. The abstract body should use 10~point type, with
a vertical spacing of 11~points, and should be indented 0.25~inches
more than normal on left-hand and right-hand margins. Insert
0.4~inches of blank space after the body. Keep your abstract brief and 
self-contained,
limiting it to one paragraph and roughly 4--6 sentences.  Gross violations will require correction at the camera-ready phase.

\subsection{Partitioning the Text} 

You should organize your paper into sections and paragraphs to help
readers place a structure on the material and understand its
contributions.

\subsubsection{Sections and Subsections}

Section headings should be numbered, flush left, and set in 11~pt bold
type with the content words capitalized. Leave 0.25~inches of space
before the heading and 0.15~inches after the heading.

Similarly, subsection headings should be numbered, flush left, and set
in 10~pt bold type with the content words capitalized. Leave
0.2~inches of space before the heading and 0.13~inches afterward.

Finally, subsubsection headings should be numbered, flush left, and
set in 10~pt small caps with the content words capitalized. Leave
0.18~inches of space before the heading and 0.1~inches after the
heading. 

Please use no more than three levels of headings.

\subsubsection{Paragraphs and Footnotes}

Within each section or subsection, you should further partition the
paper into paragraphs. Do not indent the first line of a given
paragraph, but insert a blank line between succeeding ones.
 
You can use footnotes\footnote{For the sake of readability, footnotes
should be complete sentences.} to provide readers with additional
information about a topic without interrupting the flow of the paper. 
Indicate footnotes with a number in the text where the point is most
relevant. Place the footnote in 9~point type at the bottom of the
column in which it appears. Precede the first footnote in a column
with a horizontal rule of 0.8~inches.\footnote{Multiple footnotes can
appear in each column, in the same order as they appear in the text,
but spread them across columns and pages if possible.}

\begin{figure}[ht]
\vskip 0.2in
\begin{center}
\centerline{\includegraphics[width=\columnwidth]{icml_numpapers}}
\caption{Historical locations and number of accepted papers for International
  Machine Learning Conferences (ICML 1993 -- ICML 2008) and
  International Workshops on Machine Learning (ML 1988 -- ML
  1992). At the time this figure was produced, the number of
  accepted papers for ICML 2008 was unknown and instead estimated.}
\label{icml-historical}
\end{center}
\vskip -0.2in
\end{figure} 

\subsection{Figures}
 
You may want to include figures in the paper to help readers visualize
your approach and your results. Such artwork should be centered,
legible, and separated from the text. Lines should be dark and at
least 0.5~points thick for purposes of reproduction, and text should
not appear on a gray background.

Label all distinct components of each figure. If the figure takes the
form of a graph, then give a name for each axis and include a legend
that briefly describes each curve. Do not include a title inside the
figure; instead, the caption should serve this function.

Number figures sequentially, placing the figure number and caption
{\it after\/} the graphics, with at least 0.1~inches of space before
the caption and 0.1~inches after it, as in
Figure~\ref{icml-historical}.  The figure caption should be set in
9~point type and centered unless it runs two or more lines, in which
case it should be flush left.  You may float figures to the top or
bottom of a column, and you may set wide figures across both columns
(use the environment {\tt figure*} in \LaTeX), but always place
two-column figures at the top or bottom of the page.

\subsection{Algorithms}

If you are using \LaTeX, please use the ``algorithm'' and ``algorithmic'' 
environments to format pseudocode. These require 
the corresponding stylefiles, algorithm.sty and 
algorithmic.sty, which are supplied with this package. 
Algorithm~\ref{alg:example} shows an example. 

\begin{algorithm}[tb]
   \caption{Bubble Sort}
   \label{alg:example}
\begin{algorithmic}
   \STATE {\bfseries Input:} data $x_i$, size $m$
   \REPEAT
   \STATE Initialize $noChange = true$.
   \FOR{$i=1$ {\bfseries to} $m-1$}
   \IF{$x_i > x_{i+1}$} 
   \STATE Swap $x_i$ and $x_{i+1}$
   \STATE $noChange = false$
   \ENDIF
   \ENDFOR
   \UNTIL{$noChange$ is $true$}
\end{algorithmic}
\end{algorithm}
 
\subsection{Tables} 
 
You may also want to include tables that summarize material. Like 
figures, these should be centered, legible, and numbered consecutively. 
However, place the title {\it above\/} the table with at least 
0.1~inches of space before the title and the same after it, as in 
Table~\ref{sample-table}. The table title should be set in 9~point 
type and centered unless it runs two or more lines, in which case it
should be flush left.

% Note use of \abovespace and \belowspace to get reasonable spacing 
% above and below tabular lines. 

\begin{table}[t]
\caption{Classification accuracies for naive Bayes and flexible 
Bayes on various data sets.}
\label{sample-table}
\vskip 0.15in
\begin{center}
\begin{small}
\begin{sc}
\begin{tabular}{lcccr}
\hline
\abovespace\belowspace
Data set & Naive & Flexible & Better? \\
\hline
\abovespace
Breast    & 95.9$\pm$ 0.2& 96.7$\pm$ 0.2& $\surd$ \\
Cleveland & 83.3$\pm$ 0.6& 80.0$\pm$ 0.6& $\times$\\
Glass2    & 61.9$\pm$ 1.4& 83.8$\pm$ 0.7& $\surd$ \\
Credit    & 74.8$\pm$ 0.5& 78.3$\pm$ 0.6&         \\
Horse     & 73.3$\pm$ 0.9& 69.7$\pm$ 1.0& $\times$\\
Meta      & 67.1$\pm$ 0.6& 76.5$\pm$ 0.5& $\surd$ \\
Pima      & 75.1$\pm$ 0.6& 73.9$\pm$ 0.5&         \\
\belowspace
Vehicle   & 44.9$\pm$ 0.6& 61.5$\pm$ 0.4& $\surd$ \\
\hline
\end{tabular}
\end{sc}
\end{small}
\end{center}
\vskip -0.1in
\end{table}

Tables contain textual material that can be typeset, as contrasted 
with figures, which contain graphical material that must be drawn. 
Specify the contents of each row and column in the table's topmost
row. Again, you may float tables to a column's top or bottom, and set
wide tables across both columns, but place two-column tables at the
top or bottom of the page.
 
\subsection{Citations and References} 

Please use APA reference format regardless of your formatter
or word processor. If you rely on the \LaTeX\/ bibliographic 
facility, use {\tt natbib.sty} and {\tt icml2016.bst} 
included in the style-file package to obtain this format.

Citations within the text should include the authors' last names and
year. If the authors' names are included in the sentence, place only
the year in parentheses, for example when referencing Arthur Samuel's
pioneering work \yrcite{Samuel59}. Otherwise place the entire
reference in parentheses with the authors and year separated by a
comma \cite{Samuel59}. List multiple references separated by
semicolons \cite{kearns89,Samuel59,mitchell80}. Use the `et~al.'
construct only for citations with three or more authors or after
listing all authors to a publication in an earlier reference \cite{MachineLearningI}.

Authors should cite their own work in the third person
in the initial version of their paper submitted for blind review.
Please refer to Section~\ref{author info} for detailed instructions on how to
cite your own papers.

Use an unnumbered first-level section heading for the references, and 
use a hanging indent style, with the first line of the reference flush
against the left margin and subsequent lines indented by 10 points. 
The references at the end of this document give examples for journal
articles \cite{Samuel59}, conference publications \cite{langley00}, book chapters \cite{Newell81}, books \cite{DudaHart2nd}, edited volumes \cite{MachineLearningI}, 
technical reports \cite{mitchell80}, and dissertations \cite{kearns89}. 

Alphabetize references by the surnames of the first authors, with
single author entries preceding multiple author entries. Order
references for the same authors by year of publication, with the
earliest first. Make sure that each reference includes all relevant
information (e.g., page numbers).

\subsection{Software and Data}

We strongly encourage the publication of software and data with the
camera-ready version of the paper whenever appropriate.  This can be
done by including a URL in the camera-ready copy.  However, do not
include URLs that reveal your institution or identity in your
submission for review.  Instead, provide an anonymous URL or upload
the material as ``Supplementary Material'' into the CMT reviewing
system.  Note that reviewers are not required to look a this material
when writing their review.

% Acknowledgements should only appear in the accepted version. 
\section*{Acknowledgements} 
 
\textbf{Do not} include acknowledgements in the initial version of
the paper submitted for blind review.

If a paper is accepted, the final camera-ready version can (and
probably should) include acknowledgements. In this case, please
place such acknowledgements in an unnumbered section at the
end of the paper. Typically, this will include thanks to reviewers
who gave useful comments, to colleagues who contributed to the ideas, 
and to funding agencies and corporate sponsors that provided financial 
support.  


% In the unusual situation where you want a paper to appear in the
% references without citing it in the main text, use \nocite
\nocite{langley00}

\bibliography{example_paper}
\bibliographystyle{icml2016}

\end{document} 


% This document was modified from the file originally made available by
% Pat Langley and Andrea Danyluk for ICML-2K. This version was
% created by Lise Getoor and Tobias Scheffer, it was slightly modified  
% from the 2010 version by Thorsten Joachims & Johannes Fuernkranz, 
% slightly modified from the 2009 version by Kiri Wagstaff and 
% Sam Roweis's 2008 version, which is slightly modified from 
% Prasad Tadepalli's 2007 version which is a lightly 
% changed version of the previous year's version by Andrew Moore, 
% which was in turn edited from those of Kristian Kersting and 
% Codrina Lauth. Alex Smola contributed to the algorithmic style files.  

notation.tex
\subsection{Notation}


prelim.tex
\section{Preliminaries}
\label{sec:prelim}
We begin by introducing some notations.
For any integer $k$ we use $[k]$ to denote the set of integers from $1$ to $k$.
We use lower case bold letters  for vectors (e.g.\ $\vec{v} \in \R^{d_1}$),
upper case bold letters for matrices (e.g.\ $\M \in \R^{d_1 \times d_2}$) and
bold calligraphic letters for higher order tensors (e.g.\ $\T \in \R^{d_1
\times d_2 \times d_3}$).
%
The identity matrix will be written as $\I$.
%
The $i$th row (resp. column) of a matrix $\M$ will be denoted by
$\M_{i,:}$ (resp. $\M_{:,i}$). This notation is extended to
slices of a tensor in the straightforward way.
%
If $\vec{v} \in \R^{d_1}$ and $\vec{v}' \in \R^{d_2}$, we use $\vec{v} \kron \vec{v}' \in \R^{d_1
\cdot d_2}$ to denote the Kronecker product between vectors, and its
straightforward extension to matrices and tensors.
%
Given a matrix $\M \in \R^{d_1 \times d_2}$, we use $\vectorize{\M} \in \R^{d_1
\cdot d_2}$ to denote the column vector obtained by concatenating the columns of
$\M$.

\subsection{Tensors and Tucker Decomposition}
We first recall some basic definitions of tensor algebra; more details can be found
in~\citet{Kolda09}. 
A \emph{tensor} $\T\in \R^{d_1\times\cdots \times d_p}$ is a collection
of real numbers $(\T_{i_1,\cdots,i_p}\ : \ i_n\in [d_n], n\in [p])$. The
\emph{mode-$n$} fibers of $\T$ are the vectors obtained by fixing all
indices except for the $n$th one, e.g. $\T_{:,i_2,\cdots,i_p}\in\R^{d_1}$.
The \emph{$n$th mode matricization} of $\T$ is the matrix having the
mode-$n$ fibers of $\T$ for columns and is denoted by
$\tenmat{T}{n}\in \R^{d_n\times d_1\cdots d_{n-1}d_{n+1}\cdots d_p}$.
The vectorization of a tensor is defined by $\vectorize{\T}=\vectorize{\tenmat{T}{1}}$.
The \emph{inner product} between two tensors $\ten{S}$ and $\ten{T}$ (of the same size)
is defined by $\langle \ten{S},\T\rangle = \langle\vectorize{\ten{S}},\vectorize{\T}\rangle$
and the Frobenius norm is defined by $\norm{\T}^2_F = \langle \T,\T\rangle$.
In the following $\T$ always denotes a tensor of size $d_1\times\cdots \times d_p$.

\paragraph{mode-$n$ product.} The \emph{mode-$n$ matrix product} of the tensor $\T$ and a matrix
$\X\in\R^{m\times d_n}$ is a tensor  denoted by $\T\ttm{n}\X$. It is 
of size $d_1\times\cdots \times d_{n-1}\times m \times d_{n+1}\times
\cdots \times d_p$ and is defined by the relation 
$\Y = \T\ttm{n}\X \Leftrightarrow \tenmat{Y}{n} = \X\tenmat{T}{n}$.
The \emph{mode-$n$ vector product} of the tensor $\T$ and a vector
$\vec{v}\in\R^{d_n}$ is a tensor defined by $\T\ttv{n}\vec{v} = \T\ttm{n}\vec{v}^\top
\in \R^{d_1\times\cdots \times d_{n-1}\times d_{n+1}\times
\cdots \times d_p}$.
Given tensors $\ten{S}$ and $\T$ and matrices $\X,\mat{A}$ and $\mat{B}$,
it is easy to check that
$\langle \ten{T}\ttm{n}\X, \ten{S} \rangle= \langle \ten{T},\ten{S}\ttm{n}\X^\top\rangle$ and 
$(\T\ttm{n}\mat{A})\ttm{n}\mat{B} = \T\ttm{n}\mat{BA}$ (where we assumed conforming
dimensions of the tensors and matrices).

\paragraph{Multilinear rank.} The \emph{mode-$n$ rank} of $\T$ is the dimension of the space
spanned by its mode-$n$ fibers, that is $\rank_n(\T) = \rank(\tenmat{T}{n})$. 
The \emph{multilinear rank} of $\T$, denoted by $\rank(\T)$, is the tuple of mode-$n$
 ranks of $\T$:  $rank(\T) = (R_1,\cdots,R_p)$ where $R_n = \rank_n(\T)$ for $n\in[p]$.
 We will write $\rank(\T)\leq (S_1,\cdots,S_p)$ whenever 
 $\rank_1(\T)\leq S_1,\rank_2(\T)\leq S_2,\cdots,\rank_p(\T)\leq S_p$.
 
\paragraph{Tucker decomposition.} The \emph{Tucker decomposition} is a form of
higher-order principal component analysis. It decomposes a tensor $\T$ into a
core tensor $\G$ transformed by an orthogonal matrix along each mode:
\begin{equation}
\label{eq:tucker}
\T = \G \ttm{1}\U_1\ttm{2}\U_2\ttm{3}\cdots\ttm{p}\U_p\enspace,
\end{equation}
where $\G\in \R^{R_1\times R_2 \times \cdots \times R_p}$, $\U_i \in \R^{d_i\times R_i}$ for $i\in[p]$
and $\U_i^\top \U_i = \I$ for all $i\in [p]$. The number of parameters involved in a Tucker decomposition
can be considerably smaller than $d_1d_2\cdots d_p$.
We have the following identities when matricizing and
vectorizing a Tucker decomposition
\begin{align}
\hspace{-0.1cm}\tenmat{T}{n} &= \U_n \tenmat{G}{n}(\U_p\kron\cdots\kron\U_{n+1}\kron\U_{n-1}\kron\cdots\kron\U_1)^\top\nonumber \\
&\vectorize{\T} = (\U_p\kron \U_{p-1}\kron \cdots\kron\U_1)\vectorize{\G}\enspace.\label{eq:vecTucker}
\end{align}

It is well known that $\T$ admits the Tucker 
decomposition~(\ref{eq:tucker}) if and only if  $\rank(\T)\leq  (R_1,\cdots,R_p)$ (see e.g.~\citet{Kolda09}).
Finding an exact Tucker decomposition can be done using the higher-order SVD algorithm (HOSVD) introduced
by \citet{de2000multilinear}.
Although finding the best approximation of multilinear rank $(R_1,\cdots,R_p)$ of a tensor $\T$ is a difficult
problem,  the truncated HOSVD algorithm provides good approximation guarantees 
and often performs well in practice.
	

\subsection{Low-Rank Regression}
\label{sec:LRR}
Multivariate regression is the task of recovering a  function $f: \R^d \to \R^p$ from a set of input-output 
pairs $\{(\vec{x}^{(n)}, \vec{y}^{(n)})\}_{n=1}^N$, where the outputs are sampled from the model with an additive 
noise $\vec{y} = f(\vec{x}) + \vecs{\varepsilon}$, where $\vecs{\varepsilon}$ is the error term. 
%
To solve this problem, the \emph{ordinary least squares}~(OLS) approach assumes a linear dependence between
 input and output data and boils down to finding a matrix $\vec{W} \in \R^{d \times p}$ that minimizes the 
 squared error $\| \vec{XW} - \vec{Y}\|^2_F$, where $\vec{X} \in \R^{N\times d}$ and $\vec{Y} \in \R^{N\times p}$  
 denote the input and  the output matrices. To prevent overfitting and to avoid
 numerical instabilities a ridge regularization term (i.e. $\gamma\norm{\mat{W}}^2_F$) is often added to
 the  objective function, leading to the \emph{regularized least squares}~(RLS) method.
 The main drawback of this method is that it does not take into account the dependencies between the components 
 of the response. It performs poorly when the outputs are correlated and the true dimension of the response is less than $p$.
  Indeed, it is easy to see that the OLS/RLS approach in the multivariate setting is equivalent to performing 
  $p$ linear regressions for each scalar output $\{\vec{y}_j\}_{j=1}^p$ independently.
%

\emph{Low-rank regression} (or reduced-rank regression) addresses this issue by solving the following rank penalized problem
\begin{equation}
\label{eq:matrixRRR}
 \min_{\mat{W}\in \R^{d\times p}} \| \vec{XW} - \vec{Y}\|^2_F + \gamma\norm{\mat{W}}^2_F 
  \mbox{ s.t. } \rank(\vec{W}) \leq R\enspace,
 \end{equation}
for a given integer $R$. 

Suppose that $f(\x)=\mat{W}\x$. The constraint $rank(\vec{W}) = R$ implies linear restrictions on $\vec{W}$, i.e. there must exist $p-R$ linearly 
independent vectors $\v_1, \cdots, \v_{p-R} \in \R^p$ such that $\vec{W}\v_i = \vec{0}$, 
which, in turn, implies that 
$\vec{y}^\top \v_i =  \vecs{\varepsilon}^\top \v_i$ 
for each $i\in[p-R]$. Intuitively this means that the linear subspace generated by the $\v_i$'s only contains noise. 
Another way to see this is to write $\vec{W} = \vec{AB}$ with $\vec{A}\in \R^{d\times R}$ and $\vec{B} \in \R^{R\times p}$, 
which implies that the relation between $\vec{x}$ and $\vec{y}$ is explained in an $R$-dimensional latent space. 

The rank constraint in~(\ref{eq:matrixRRR})  was first proposed in~\citep{Anderson51}, 
whereas the term \emph{reduced-rank regression} was introduced in~\citep{Izenman75}.
Adding a  ridge regularization to the rank penalized problem was proposed in~\cite{Mukherjee2011}. 
In the rest of the paper we will refer to this approach as low-rank regression (LRR).
The solution of  minimization
problem~(\ref{eq:matrixRRR}) is given by projecting the RLS solution onto the space spanned by the top $R$
eigenvectors of $\mat{Y}^\top\mat{P}\mat{Y}$ where $\mat{P}$ is the orthogonal projection 
matrix onto the column space of $\X$, that is $\mat{W}_{LRR} = \mat{W}_{RLS}\mats{\Pi}$ where
$\mats{\Pi}$ is the matrix of the aforementioned projection.  
For more description and discussion of reduced-rank regression, we refer
the reader to the books of~\citet{Reinsel1998} and~\citet{Izenman08}.

%
%The rank constraint in RRR can be interpreted from several perspectives. 
%First, the constraint $rank(\vec{W}) = R$ implies linear restrictions on $\vec{W}$, i.e. there must exist $p-R$ linearly 
%independent vectors $\vecs{\ell}_1, \cdots, \vecs{\ell}_{p-R} \in \R^p$ such that $\vec{W}\vecs{\ell}_i = \vec{0}$, 
%which implies that 
%$\vec{y}^\top \vecs{\ell}_i =  (\vec{x}^\top \vec{W} + \vecs{\varepsilon}^\top) \vecs{\ell}_i =  \vecs{\varepsilon}^\top \vecs{\ell}_i$ 
%for each $i$. Second, we can write $\vec{W} = \vec{AB}$ with $\vec{A}\in \R^{d\times R}$ and $\vec{B} \in \R^{R\times p}$, 
%which implies that the relation between $\vec{x}$ and $\vec{y}$ is explained in an $R$-dimensional latent space. 
%Third, if we consider the linear regression function $f_j:\R^d\to\R$, for $j\in \{1,\ldots,p\}$, 
%that predicts the $j$th component of the response $\vec{y}\in \mathbb{R}^p$ (i.e. $f_j(\vec{x}) = \sum_{i=1} \vec{x}_i\vec{W}_{ij}$), 
%the rank constraint implies that the dimension of the functional space spanned by the $f_j$'s is less than $R$. 




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%






proofs.tex
Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. Here goes the proofs. 
rademacher.tex


We now provide a generalization bound of our algorithm based on the notion of the rademacher complexity.


\begin{definition}
For $n\in\{ 1,\ldots, N\}$ and $m_i\in\{ 1,\ldots, d_i\}$, $1\leq i \leq p$, let $\sigma^n_{m_1 m_2 \ldots m_p}$ be independent random variables, distributed uniformly in $\{ -1, 1\}$. The empirical Rademacher complexity
of a class $\mathcal{F}$ of tensor-valued functions ${f}:\mathbb{R}^{d_0} \to \mathbb{R}^{d_1\times d_2 \times \ldots \times d_p}$ is the function $\hat{\mathcal{R}}(\mathcal{F})$ defined on $\mat{x} = (x_n)_{n=1}^N \in \mathbb{R}^{d_0N}$ by
\begin{equation*}
\hat{\mathcal{R}}(\mathcal{F})(\mat{x}) = \frac{1}{N} \mathbb{E}_\ten{\sigma} \left[ \sup_{{f}\in\mathcal{F}}  \sum_{n=1}^N \langle \ten{\sigma}^{(n)} , {f}(x_n) \rangle   \right],
\end{equation*}
where $\ten{\sigma}$ is the tensor $\big(\sigma^n_{m_1 m_2 \ldots m_p}\big)\in \mathbb{R}^{N \times d_1 \times \ldots \times d_p}$.
\end{definition}



\begin{theorem}
Let $\mathcal{F}$ be a set of linear tensor functions with Tucker rank penalization, i.e., $\mathcal{F} = \left\{ x\mapsto \W\ttm 1 x : rank(\W) = (R_0,R_1,\ldots,R_p) \right\}$ with $\W\in\mathbb{R}^{d_0\times d_1 \times \ldots \times d_p}$ and $\forall x\in\mathbb{R}^{d_0}$, and let $\mat{x}\in\mathbb{R}^{d_0N}$ with $\|x_n\|^2\leq 1$.
%
The empirical Rademacher complexity of $\mathcal{F}$ can be upper bounded as follows:
\begin{equation*}
\hat{\mathcal{R}}(\mathcal{F})(\mat{x}) \leq \sqrt{\frac{R_0\times R_1\times\ldots\times R_p}{N}} \sup_{\G}\|\G\|_{max},
\end{equation*}
wehre $\G\in\mathbb{R}^{R_0\times\ldots\times R_p}$ is the core tensor of $\W$.

\end{theorem}

\begin{proof}
\begin{align*}
\hat{\mathcal{R}}(\mathcal{F})(\mat{x}) &= \frac{1}{N} \mathbb{E}_\ten{\sigma} \left[ \sup_{{f}\in\mathcal{F}}  \sum_{n=1}^N \langle \ten{\sigma}^{(n)} , {f}(x_n) \rangle   \right]\\
%
&=  \frac{1}{N} \mathbb{E}_\ten{\sigma} \left[ \sup_{\W\in\mathbb{R}^{d_0\times  \ldots \times d_p}}  \sum_{n=1}^N \langle \ten{\sigma}^{(n)} , \W\ttm 1x_n \rangle   \right]\\
&=  \frac{1}{N} \mathbb{E}_\ten{\sigma} \left[ \sup_{\W}  \sum_{n=1}^N \langle x_n \otimes \ten{\sigma}^{(n)} , \W \rangle   \right]\\
%
& \qquad (Cauchy-Schwarz)\\
&\leq  \frac{1}{N}  \sup_{\W} \|\W\| \text{\ }  \mathbb{E}_\ten{\sigma} \left[ \Big\| \sum_{n=1}^N  x_n \otimes \ten{\sigma}^{(n)}\Big\|   \right]\\
%
& \qquad (Jensen's \ inequality)\\
&\leq  \frac{1}{N}  \sup_{\W} \|\W\| \text{\ }  \left[ \mathbb{E}_\ten{\sigma} \bigg[ \Big\| \sum_{n=1}^N  x_n \otimes \ten{\sigma}^{(n)}\Big\|^2 \bigg]  \right]^{\frac{1}{2}}\\
%
& \qquad (?)\\
&=  \frac{1}{N}  \sup_{\W} \|\W\| \text{\ }  \Bigg[ \mathbb{E}_\ten{\sigma} \bigg[ \Big\| \mat{X}^{\top} \ten{\sigma}_{(1)}\Big\|^2 \bigg]  \Bigg]^{\frac{1}{2}}\\
%
&=  \frac{1}{N}  \sup_{\W} \|\W\| \text{\ }  \Bigg[ \mathbb{E}_\ten{\sigma} \bigg[ tr\Big(\ten{\sigma}_{(1)} ^{\top}  \mat{X} \mat{X}^{\top}\ten{\sigma}_{(1)}\Big)  \bigg]  \Bigg]^{\frac{1}{2}}\\
%
&=  \frac{1}{N}  \sup_{\W} \|\W\| \text{\ }  \Bigg[  tr\bigg(  \mathbb{E}_\ten{\sigma} \Big[\ten{\sigma}_{(1)}\ten{\sigma}_{(1)}^{\top}\Big]  \mat{X} \mat{X}^{\top}\bigg)    \Bigg]^{\frac{1}{2}}\\
%
&\qquad \big(\mathbb{E}_\ten{\sigma} [\ten{\sigma}_{(1)}\ten{\sigma}_{(1)}^{\top}]   = I\big) \text{\ by the independence and unit-variance of the Rademacher variables}\\
&=  \frac{1}{N}  \sup_{\W} \|\W\| \text{\ }  \bigg[ tr\Big(\mat{X} \mat{X}^{\top}\Big)  \bigg] ^{\frac{1}{2}}\\
%
& \qquad (\forall 1\leq n\leq N, \|x_n\|^2\leq 1 \Rightarrow tr\Big(\mat{X} \mat{X}^{\top}\Big) \leq N)\\
&\leq  \frac{1}{\sqrt{N}}  \sup_{\W} \|\W\| \text{\ }\\
& \qquad (\| W \| = \| \G \|) \text{\textcolor{blue}{ajouter U orthogonaux et decomposition tucker}}\\
&=  \frac{1}{\sqrt{N}}  \sup_{\G\in\mathbb{R}^{R_0\times\ldots\times R_p}} \|\G\| \text{\ }\\
%
&\leq \sqrt{\frac{R_0\times R_1\times\ldots\times R_p}{N}} \sup_{\G}\|\G\|_{max}\\
\end{align*}

\textcolor{blue}{verifier notation: $\mat{X}\in\mathbb{R}^{N\times d_0}$ dans la preuve et $\mat{x}\in\mathbb{R}^{d_0 N}$ dans le theorem et definition}

\end{proof}



supmat.tex
\begin{theorem}
Let $\mathcal{F}$ be a set of linear tensor functions with Tucker rank penalization, i.e., $\mathcal{F} = \left\{ x\mapsto \W\ttm 1 x : rank(\W) = (R_0,R_1,\ldots,R_p) \right\}$ with $\W\in\mathbb{R}^{d_0\times d_1 \times \ldots \times d_p}$ and $\forall x\in\mathbb{R}^{d_0}$, and let $\mat{x}\in\mathbb{R}^{d_0N}$ with $\|x_n\|^2\leq 1$.
%
The empirical Rademacher complexity of $\mathcal{F}$ can be upper bounded as follows:
\begin{equation*}
\hat{\mathcal{R}}(\mathcal{F})(\mat{x}) \leq \sqrt{\frac{R_0\times R_1\times\ldots\times R_p}{N}} \sup_{\G}\|\G\|_{max},
\end{equation*}
wehre $\G\in\mathbb{R}^{R_0\times\ldots\times R_p}$ is the core tensor of $\W$.

\end{theorem}

\begin{proof}
\begin{align*}
\hat{\mathcal{R}}(\mathcal{F})(\mat{x}) &= \frac{1}{N} \mathbb{E}_\ten{\sigma} \left[ \sup_{{f}\in\mathcal{F}}  \sum_{n=1}^N \langle \ten{\sigma}^{(n)} , {f}(x_n) \rangle   \right]\\
%
&=  \frac{1}{N} \mathbb{E}_\ten{\sigma} \left[ \sup_{\W\in\mathbb{R}^{d_0\times  \ldots \times d_p}}  \sum_{n=1}^N \langle \ten{\sigma}^{(n)} , \W\ttm 1x_n \rangle   \right]\\
&=  \frac{1}{N} \mathbb{E}_\ten{\sigma} \left[ \sup_{\W}  \sum_{n=1}^N \langle x_n \otimes \ten{\sigma}^{(n)} , \W \rangle   \right]\\
%
& \qquad (Cauchy-Schwarz)\\
&\leq  \frac{1}{N}  \sup_{\W} \|\W\| \text{\ }  \mathbb{E}_\ten{\sigma} \left[ \Big\| \sum_{n=1}^N  x_n \otimes \ten{\sigma}^{(n)}\Big\|   \right]\\
%
& \qquad (Jensen's \ inequality)\\
&\leq  \frac{1}{N}  \sup_{\W} \|\W\| \text{\ }  \left[ \mathbb{E}_\ten{\sigma} \bigg[ \Big\| \sum_{n=1}^N  x_n \otimes \ten{\sigma}^{(n)}\Big\|^2 \bigg]  \right]^{\frac{1}{2}}\\
%
& \qquad (?)\\
&=  \frac{1}{N}  \sup_{\W} \|\W\| \text{\ }  \Bigg[ \mathbb{E}_\ten{\sigma} \bigg[ \Big\| \mat{X}^{\top} \ten{\sigma}_{(1)}\Big\|^2 \bigg]  \Bigg]^{\frac{1}{2}}\\
%
&=  \frac{1}{N}  \sup_{\W} \|\W\| \text{\ }  \Bigg[ \mathbb{E}_\ten{\sigma} \bigg[ tr\Big(\ten{\sigma}_{(1)} ^{\top}  \mat{X} \mat{X}^{\top}\ten{\sigma}_{(1)}\Big)  \bigg]  \Bigg]^{\frac{1}{2}}\\
%
&=  \frac{1}{N}  \sup_{\W} \|\W\| \text{\ }  \Bigg[  tr\bigg(  \mathbb{E}_\ten{\sigma} \Big[\ten{\sigma}_{(1)}\ten{\sigma}_{(1)}^{\top}\Big]  \mat{X} \mat{X}^{\top}\bigg)    \Bigg]^{\frac{1}{2}}\\
%
&\qquad \big(\mathbb{E}_\ten{\sigma} [\ten{\sigma}_{(1)}\ten{\sigma}_{(1)}^{\top}]   = I\big) \text{\ by the independence and unit-variance of the Rademacher variables}\\
&=  \frac{1}{N}  \sup_{\W} \|\W\| \text{\ }  \bigg[ tr\Big(\mat{X} \mat{X}^{\top}\Big)  \bigg] ^{\frac{1}{2}}\\
%
& \qquad (\forall 1\leq n\leq N, \|x_n\|^2\leq 1 \Rightarrow tr\Big(\mat{X} \mat{X}^{\top}\Big) \leq N)\\
&\leq  \frac{1}{\sqrt{N}}  \sup_{\W} \|\W\| \text{\ }\\
& \qquad (\| W \| = \| \G \|) \text{\textcolor{blue}{ajouter U orthogonaux et decomposition tucker}}\\
&=  \frac{1}{\sqrt{N}}  \sup_{\G\in\mathbb{R}^{R_0\times\ldots\times R_p}} \|\G\| \text{\ }\\
%
&\leq \sqrt{\frac{R_0\times R_1\times\ldots\times R_p}{N}} \sup_{\G}\|\G\|_{max}\\
\end{align*}

\textcolor{blue}{verifier notation: $\mat{X}\in\mathbb{R}^{N\times d_0}$ dans la preuve et $\mat{x}\in\mathbb{R}^{d_0 N}$ dans le theorem et definition}

\end{proof}



supplementary.tex
\documentclass{article}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ICML default packages
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% use Times
\usepackage{times}
% For figures
\usepackage{graphicx} % more modern
%\usepackage{epsfig} % less modern
\usepackage{subfigure} 

% For citations
\usepackage{natbib}

% For algorithms
\usepackage{algorithm}
\usepackage{algorithmic}

\usepackage{hyperref}

\newcommand{\theHalgorithm}{\arabic{algorithm}}

\usepackage{icml2016} 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% our packages
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{mathtools}
\usepackage{lmodern}

% for notes (example: \hachem{We need to change all the notations...} \guillaume{No way!})
\usepackage[textsize=footnotesize,color=green!40]{todonotes}

% for \llbracket and \rrbracket
\usepackage{stmaryrd}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% our macros
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\input{macros.tex}

% The \icmltitle you define below is probably too long as a header.
% Therefore, a short form for the running title is supplied here:
\icmltitlerunning{Higher-Order Low-Rank Regression (Supplementary Material)}


\begin{document}
\twocolumn[
\icmltitle{Higher-Order Low-Rank Regression (Supplementary Material)}

% It is OKAY to include author information, even for blind
% submissions: the style file will automatically remove it for you
% unless you've provided the [accepted] option to the icml2016
% package.
\icmlauthor{Your Name}{email@yourdomain.edu}
\icmladdress{Your Fantastic Institute,
            314159 Pi St., Palo Alto, CA 94306 USA}
\icmlauthor{Your CoAuthor's Name}{email@coauthordomain.edu}
\icmladdress{Their Fantastic Institute,
            27182 Exp St., Toronto, ON M6H 2T1 CANADA}

% You may provide any keywords that you 
% find helpful for describing your paper; these are used to populate 
% the "keywords" metadata in the PDF but will not be shown in the document
\icmlkeywords{boring formatting information, machine learning, ICML}

\vskip 0.3in

\input{proofs}

\bibliographystyle{icml2016}
\bibliography{bib_TensorKernels}

]
\end{document}



xp.tex

 \section{Experiments}
\label{sec:xp}
We present and analyze the experimental results obtained on
a regression task on synthetic data, an
image reconstruction task and a meteorological forecasting task on real data%
\footnote{The code and data used in the experiments will be made available by the authors.}.
We compare the predictive accuracy of  HOLRR with the following methods:\newline
- RLS: Regularized least squares.\newline
- LRR: Low-rank regression (see Section~\ref{sec:LRR}).\newline
- ADMM: a multilinear approach based on tensor trace norm regularization introduced in \citep{gandy2011tensor}
 and \citep{romera2013}.\newline
 - MLMT-NC: a nonconvex approach proposed in \citep{romera2013}  in the context of multilinear multitask learning.
 
 For experiments with kernel algorithms we use the readily available kernelized RLS
  and the LRR kernel extension  proposed 
 in~\citep{Mukherjee2011} (note that ADMM and MLMT-NC only consider a linear
 dependency between inputs and outputs).
 
 Although  MLMLT-NC is perhaps the closest algorithm to ours, we applied it only to simulated data. 
 This is because MLMLT-NC is computationally  very expensive and becomes intractable for large data sets.
  To give an idea on the running times,
 averages for some of the experiments are reported in Table~\ref{tab:runtime}.
 
 \begin{table}
 \begin{center}
 
 \caption{Average running time for different size of the regression tensor $\W$. Each line
 comes from one of the experiments. The times reported
 are the average learning times in seconds for a training set of size $100$.}
 \label{tab:runtime}
  \begin{scriptsize}
  
 \vspace{0.1cm}
 
 \begin{sc}
 \setlength{\tabcolsep}{3pt}
 \begin{tabular}{c|ccccc}
 \hline
size of $\W$ & RLS & LRR & ADMM & MLMT-NC & HOLRR\\

\hline
\rule{0pt}{2.6ex}

${10\times 10\times 10\times 10}$ &$0.012$ & $0.45$ & $12.92$ & $945.79 $ & $0.04$ \\
${3\times 70\times 70}$ & $0.04$ & $10.46$ & $-$&$-$&$0.08$\\ 
${160\times 5\times 16 \times 5}$  & $2.07$ & $2.24$ & $40.23$ & $-$ & $1.67$ \\
\hline
 \end{tabular}
 \end{sc}
 \end{scriptsize}
 
 \vspace{-0.5cm}
 
 
 \end{center}
 \end{table}
 
\subsection{Synthetic Data}
We generate both linear and nonlinear data.
The linear data is drawn from the model  $\ten{Y} = \ten{W} \ttv{1} \vec{x} + \ten{E}$
where $\ten{W}\in \R^{10\times 10\times 10 \times 10}$ is a tensor of multilinear 
rank $(6,4,4,8)$ drawn at random, $\x\in\R^{10}$ is drawn from $\mathcal{N}(0,\I)$ and each component
of the error tensor $\ten{E}$ is drawn from $\mathcal{N}(0,0.1)$. The nonlinear data
is drawn from the model $\ten{Y} = \ten{W} \ttv{1} (\vec{x}\kron\x) + \ten{E}$ where
$\ten{W}\in \R^{25\times 10\times 10 \times 10}$ is a tensor of  
rank $(5,6,4,2)$ drawn at random and $\x\in\R^5$ and $\ten{E}$ are generated as above.
The hyper-parameters for all algorithms are selected using $3$-fold cross validation
on the training data. 

These experiments have been carried out for different sizes of the training data set, 20 trials have been executed for each size.  
The average RMSEs on a test set of size 100 for the 20 trials are reported in Figure~\ref{fig:xp_synth}. We see
that the HOLRR algorithm clearly outperforms the other methods on the linear data. MLMT-NC is the method
obtaining the better accuracy after ours, it is however much more computationally expensive%
~(see Table~\ref{tab:runtime}). On the nonlinear data the LRR method achieves good performances
but HOLRR is still significantly more accurate for small data set sizes.

\begin{figure}

\begin{center}
\includegraphics[scale=0.45]{figures/XP_synth_linear.pdf}%
\includegraphics[scale=0.45,trim=1cm 0cm 0cm 0cm]{figures/XP_synth_poly_wo_ADMM.pdf}
\caption{Average RMSE as a function of the training set size: (left) linear data, (right) nonlinear data.}
\label{fig:xp_synth}
\end{center}

\vspace{-0.7cm}

\end{figure}

\subsection{Image Reconstruction from Noisy Measurements}
\label{sec:xp_images}
To give an illustrative intuition on the differences between matrix and multilinear rank regularization
 we generate data  from the model $\ten{Y} = \ten{W} \ttv{1} \vec{x} + \ten{E}$ where
the tensor $\W$ is  a color image of size $m\times n$ encoded with three color channels
 RGB. We consider two different tasks depending on the input dimension: 
 (i) $\W\in\R^{3\times m \times n}$, $\x\in \R^3$ and 
 (ii) $\W\in\R^{ n \times m \times 3}$, $\x\in \R^n$.
 In both tasks the components of both $\x$  and $\ten{E}$ are drawn from $\mathcal{N}(0,1)$ and the 
 regression tensor $\W$ is learned from a training set of size $200$. 

 This experiment allows us to visualize the tensors returned by the RLS,  LRR and HOLRR algorithms.
 The results are shown in Figure~\ref{fig:xp_images} for three images: a green cross (of size
$50\times 50$), a thumbnail of a Rothko painting ($44\times 70$) and a square made of
triangles ($70\times 70$), note that the first two images have a low rank structure which
is not the case for the third one. 

We first see that HOLRR clearly outperforms LRR on the task where the input dimension is small (task (i)).
This is to be  expected since the rank of the matrix $\tenmat{W}{1}$ is at most $3$ and
LRR is unable to enforce a low-rank structure on the output modes of $\W$. When the rank constraint is
set to $1$ for LRR and $(3,1,1)$ for HOLRR we clearly see that unlike HOLRR 
the LRR approach does not enforce any low-rank
structure on the regression tensor along the output modes.
%
On task (ii) the difference is more subtle, but we can see that setting a rank constraint
of $2$ for the LRR algorithm prevents the model from capturing the white border around the green
cross and creates the vertical lines artifact in the Rothko painting. For higher values of the rank the
model starts to learn the noise. The tensor returned by HOLRR with rank $(4,4,3)$ does not exhibit these
behaviors and gives better results
on these two images. On the square image which does not have a low-rank structure both 
algorithms do not perform very well.
%
Overall, we see that capturing the mutlilinear low-rank structure of the output data allows HOLRR to
separate the noise from the true signal better than RLS and LRR.

\begin{figure*}[th]
\vspace{-1cm}
    \hspace{-3cm}\begin{minipage}{.5\textwidth}
        \centering
        \includegraphics[height=3.5cm,,trim=0cm 0cm 0cm 0cm]{figures/XP_cross_small_din.pdf}
		\includegraphics[height=4.5cm,,trim=0cm 0.5cm 0cm 0.5cm]{figures/XP_rothko_small_din.pdf}
		\includegraphics[height=3.5cm,,trim=0cm 0cm 0cm 0cm]{figures/XP_square_small_din.pdf}
    \end{minipage}\hspace{2.5cm}%
\begin{minipage}{.5\textwidth}
        \centering
        \includegraphics[height=3.5cm,trim=0cm 0cm 0cm 0cm]{figures/XP_cross_big_din.pdf}
        %\makebox[\textwidth]{\includegraphics[=\paperwidth]{figures/XP_cross_big_din.pdf}
		\includegraphics[height=4.5cm,,trim=0cm 0.5cm 0cm 0.5cm]{figures/XP_rothko_big_din.pdf}
		\includegraphics[height=3.5cm,,trim=0cm 0cm 0cm 0cm]{figures/XP_square_big_din.pdf}
    \end{minipage}
    \caption{Image reconstruction from noisy measurements: $\Y = \W\ttv{1}\vec{x} + \ten{E}$ where $\W$ is a color image (RGB). 
    (left) Task (i): input dimension is the number of channels. (right) Task (ii): input dimension is the height
    of the image. Each image is labeled with the name of the algorithm followed by the value used for the rank constraint.}
    \label{fig:xp_images}
\end{figure*}



\subsection{Real Data}
We compare our algorithm with other methods on the task of meteo forecasting. 
We collected data from the meteorological office of the UK
\footnote{\tiny{\url{http://www.metoffice.gov.uk/public/weather/climate-historic/}}}: 
monthly average measurements of 5 variables in 16 stations across the UK from 1960 to 2000.
 The forecasting task consists in predicting the values of the 5 variables in  the 16 stations from their values in the preceding months.  
 We use the values of all the variables from the past 2 months as covariates and consider
 the tasks of predicting the values of all variables for the next $k$ months,
 the output tensors are thus in $\R^{k\times 16 \times 5}$.
 
 We randomly split the available data into a training set of size $N$, a test set of size $50$ and
 a validation set of size $20$, for all methods hyper-parameters are chosen w.r.t. to their
 performance on the validation set. The average test RMSE over $10$ runs of this experiment 
 for different values of $N$ and $k$ are reported in Table~\ref{tab:xp_meteo}. We see that the
 HOLRR has overall a better predictive accuracy than the other methods, especially when the
 tensor structure of the output data gets richer ($k=3$ and $k=5$). On this dataset, nonlinear
 methods only improve the results by a small margin with the RBF kernel, and the polynomial
 kernel performs poorly.


 
 \begin{table}
 \begin{center}
 \caption{Forecasting experiment: the task is to predict the values for all variables for the next
 $k$ months using a training set of size $N$.}
 \label{tab:xp_meteo}
 \begin{scriptsize}
 
 \vspace{0.1cm}
 
 \begin{sc}
\begin{tabular}{c|cccc}
\hline
\rule{0pt}{2.6ex}
 $k=1,N=50$ & ADMM & RLS & LRR & HOLRR \\
\hline
 linear & $0.6323$ & $0.6481$ & $\underline{0.6261}$ & $0.6310$ \\
 rbf &  - &0.6081 & 0.6076 & \textbf{\underline{0.6035}} \\
 poly &  - &0.6580 & \underline{0.6526} & 0.6793\vspace{0.1cm}\\
 \hline
\rule{0pt}{2.6ex}
$k=1,N=100$ & ADMM & RLS & LRR & HOLRR \\
\hline
linear & \underline{0.5997} & 0.6218 & \underline{0.6003} & \underline{0.6003} \\
rbf &  - &0.5926 & 0.5920 & \textbf{\underline{0.5808}} \\
poly &  - &0.6230 & \underline{0.6160} & 0.6222\vspace{0.1cm}\\
\hline
\rule{0pt}{2.6ex}
 $k=3,n=50$ & ADMM & RLS & LRR & HOLRR \\
 \hline
 linear & 0.6367 & 0.6683 & 0.6256 & \underline{0.6163} \\
 rbf &  - &0.6297 & 0.6296 & \textbf{\underline{0.6096}} \\
 poly &  - &0.6587 & 0.6412 & \underline{0.6355}\vspace{0.1cm}\\
\hline
\rule{0pt}{2.6ex}
  $k=3,N=100$ & ADMM & RLS & LRR & HOLRR \\
  \hline
  linear & 0.6269 & 0.6740 & 0.6186 & \underline{0.6086} \\
  rbf &  - &0.6153 & 0.6090 & \textbf{\underline{0.6015}} \\
  poly &  - &0.6588 & 0.6439 & \underline{0.6379}\vspace{0.1cm}\\
\hline
\rule{0pt}{2.6ex}  
$k=5,N=50$   & ADMM & RLS & LRR & HOLRR \\
\hline
linear & 0.6302 & 0.6719 & 0.6232 & \underline{0.6127} \\
rbf &  - &0.6150 & 0.6192 & \textbf{\underline{0.6017}} \\
poly &  - &0.6579 & \underline{0.6338} & 0.6353\vspace{0.1cm}\\
\hline
\rule{0pt}{2.6ex}
 $k=5,N=100$ & ADMM & RLS & LRR & HOLRR \\
 \hline
 linear & 0.6140 & 0.6449 & 0.6021 & \underline{0.5971} \\
 rbf &  - &0.6020 & 0.5908 & \textbf{\underline{0.5886}} \\
 poly &  - &0.6432 & 0.6202 & \underline{0.6107}\\
 
\hline
\end{tabular}
 
  \end{sc}
 \end{scriptsize}
 
 \end{center}
 
 
\vspace{-0.5cm}

  \end{table}
  
  
main.bbl
\begin{thebibliography}{14}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\expandafter\ifx\csname urlstyle\endcsname\relax
  \providecommand{\doi}[1]{doi: #1}\else
  \providecommand{\doi}{doi: \begingroup \urlstyle{rm}\Url}\fi

\bibitem[Anderson(1951)]{Anderson51}
Anderson, T.~W.
\newblock Estimating linear restrictions on regression coefficients for
  multivariate normal distributions.
\newblock \emph{Annals of Mathematical Statistics}, 22:\penalty0 327--351,
  1951.

\bibitem[Bahadori et~al.(2014)Bahadori, Yu, and Liu]{Bahadori2014}
Bahadori, M.~T., Yu, Q.~R., and Liu, Y.
\newblock Fast multivariate spatio-temporal analysis via low rank tensor
  learning.
\newblock In \emph{NIPS}. 2014.

\bibitem[Cichocki et~al.(2009)Cichocki, Zdunek, Phan, and Amari]{Cichocki2009}
Cichocki, A., Zdunek, R., Phan, A.H., and Amari, S.I.
\newblock \emph{{Nonnegative Matrix and Tensor Factorizations}}.
\newblock Applications to Exploratory Multi-way Data Analysis and Blind Source
  Separation. Wiley, 2009.

\bibitem[De~Lathauwer et~al.(2000)De~Lathauwer, De~Moor, and
  Vandewalle]{de2000multilinear}
De~Lathauwer, Lieven, De~Moor, Bart, and Vandewalle, Joos.
\newblock A multilinear singular value decomposition.
\newblock \emph{SIAM journal on Matrix Analysis and Applications}, 21\penalty0
  (4):\penalty0 1253--1278, 2000.

\bibitem[Gandy et~al.(2011)Gandy, Recht, and Yamada]{gandy2011tensor}
Gandy, Silvia, Recht, Benjamin, and Yamada, Isao.
\newblock Tensor completion and low-n-rank tensor recovery via convex
  optimization.
\newblock \emph{Inverse Problems}, 27\penalty0 (2):\penalty0 025010, 2011.

\bibitem[Izenman(1975)]{Izenman75}
Izenman, A.~J.
\newblock Reduced-rank regression for the multivariate linear model.
\newblock \emph{Journal of Multivariate Analysis}, 5\penalty0 (2):\penalty0
  248--264, 1975.

\bibitem[Izenman(2008)]{Izenman08}
Izenman, A.~J.
\newblock \emph{{Modern Multivariate Statistical Techniques: Regression,
  Classification, and Manifold Learning}}.
\newblock Springer-Verlag, New York, 2008.

\bibitem[Kolda \& Bader(2009)Kolda and Bader]{Kolda09}
Kolda, T.~G. and Bader, B.~W.
\newblock Tensor decompositions and applications.
\newblock \emph{SIAM review}, 51\penalty0 (3):\penalty0 455--500, 2009.

\bibitem[Lu et~al.(2013)Lu, Plataniotis, and Venetsanopoulos]{Lu2013}
Lu, H., Plataniotis, K.N., and Venetsanopoulos, A.
\newblock \emph{Multilinear Subspace Learning: Dimensionality Reduction of
  Multidimensional Data}.
\newblock CRC Press, 2013.

\bibitem[Mukherjee \& Zhu(2011)Mukherjee and Zhu]{Mukherjee2011}
Mukherjee, A. and Zhu, J.
\newblock Reduced rank ridge regression and its kernel extensions.
\newblock \emph{Statistical analysis and data mining}, 4\penalty0 (6):\penalty0
  612--622, 2011.

\bibitem[Reinsel \& Velu(1998)Reinsel and Velu]{Reinsel1998}
Reinsel, G.C. and Velu, R.P.
\newblock \emph{Multivariate reduced-rank regression: theory and applications}.
\newblock Lecture Notes in Statistics. Springer, 1998.

\bibitem[Romera-Paredes et~al.(2013)Romera-Paredes, Aung, Bianchi-Berthouze,
  and Pontil]{romera2013}
Romera-Paredes, B., Aung, M.~H., Bianchi-Berthouze, N., and Pontil, M.
\newblock Multilinear multitask learning.
\newblock In \emph{Proceedings of the 30th International Conference on Machine
  Learning (ICML-13)}, pp.\  1444--1452, 2013.

\bibitem[Signoretto et~al.(2013)Signoretto, Dinh, De~Lathauwer, and
  Suykens]{signoretto2013a}
Signoretto, M., Dinh, Q.~T., De~Lathauwer, L., and Suykens, J.~K.
\newblock Learning with tensors: a framework based on convex optimization and
  spectral regularization.
\newblock \emph{Machine Learning}, pp.\  1--49, 2013.

\bibitem[Zhou et~al.(2013)Zhou, Li, and Zhu]{Zhou2013}
Zhou, H., Li, L., and Zhu, H.
\newblock Tensor regression with applications in neuroimaging data analysis.
\newblock \emph{Journal of the American Statistical Association}, 108\penalty0
  (502):\penalty0 540--552, 2013.

\end{thebibliography}
