Alternating DCA for reduced-rank multitask linear regression with covariance matrix estimation

Abstract

We study a challenging problem in machine learning that is the reduced-rank multitask linear regression with covariance matrix estimation. The objective is to build a linear relationship between multiple output variables and input variables of a multitask learning process, taking into account the general covariance structure for the errors of the regression model in one hand, and reduced-rank regression model in another hand. The problem is formulated as minimizing a nonconvex function in two joint matrix variables (X,Θ) under the low-rank constraint on X and positive definiteness constraint on Θ. It has a double difficulty due to the non-convexity of the objective function as well as the low-rank constraint. We investigate a nonconvex, nonsmooth optimization approach based on DC (Difference of Convex functions) programming and DCA (DC Algorithm) for this hard problem. A penalty reformulation is considered which takes the form of a partial DC program. An alternating DCA and its inexact version are developed, both algorithms converge to a weak critical point of the considered problem. Numerical experiments are performed on several synthetic and benchmark real multitask linear regression datasets. The numerical results show the performance of the proposed algorithms and their superiority compared with three classical alternating/joint methods.

Appendices

Appendix A: estimate μ _Θ

The next lemma indicates how to estimate the value of μ_Θ to guarantee the convexity of function $\widetilde {H}(\cdot ,\boldsymbol {\varTheta })$ defined in Theorem 2. In this lemma, ∥⋅∥₂ denotes the ℓ₂-norm (or spectral norm) of a matrix.

Lemma 1

For each fixed $\boldsymbol {\varTheta } \in \mathbb {R}^{m \times m}$, if $\mu _{\boldsymbol {\varTheta }} \leq 1/\left (\|\frac {2}{n} {\sum }_{i=1}^{n}\boldsymbol {\phi }_{i} \boldsymbol {\phi }_{i}^{\top }\|_{2} \|\boldsymbol {\varTheta }\|_{2} + 2\alpha \right )$, then $\widetilde {H}(\cdot ,\boldsymbol {\varTheta })$ is convex.

Proof

Knowing that the function $\max \limits _{\boldsymbol {Y} \in \mathcal {X}} (2\langle \boldsymbol {X}, \boldsymbol {Y} \rangle - \|\boldsymbol {Y}\|_{F}^{2})$ is convex in X and the sum of two convex functions is also convex, then it is sufficient to choose μ_Θ such that the function $(\frac {1}{2\mu _{\boldsymbol {\varTheta }}}-\alpha ) \|\boldsymbol {X}\|_{F}^{2} - \frac {1}{n} {\sum }_{i=1}^{n} (\boldsymbol {z}_{i}-\boldsymbol {X} \boldsymbol {\phi }_{i})^{\top } \boldsymbol {\varTheta } (\boldsymbol {z}_{i}-\boldsymbol {X} \boldsymbol {\phi }_{i})$ becomes convex. For this aim, we can take $(\frac {1}{\mu _{\boldsymbol {\varTheta }}}-2\alpha )$ greater than or equal to the spectral radius of the Hessian matrix of ${\varLambda }(\boldsymbol {X}) = \frac {1}{n} {\sum }_{i=1}^{n} (\boldsymbol {z}_{i}-\boldsymbol {X} \boldsymbol {\phi }_{i})^{\top } \boldsymbol {\varTheta } (\boldsymbol {z}_{i}-\boldsymbol {X} \boldsymbol {\phi }_{i})$, i.e. $\frac {1}{\mu _{\boldsymbol {\varTheta }}}-2\alpha \geq \rho (\nabla ^{2} {\varLambda }(\boldsymbol {X}))$. From the matrix differential calculus (see, e.g., [38]), we have

$$ \begin{array}{@{}rcl@{}} \nabla {\varLambda}(\boldsymbol{X}) &=& \frac{2}{n} \sum\limits_{i=1}^{n} \boldsymbol{\varTheta}(\boldsymbol{X}\boldsymbol{\phi}_{i}-\boldsymbol{z}_{i})\boldsymbol{\phi}_{i}^{\top} = \boldsymbol{\varTheta} \boldsymbol{X} \left( \frac{2}{n} \sum\limits_{i=1}^{n}\boldsymbol{\phi}_{i}\boldsymbol{\phi}_{i}^{\top}\right) - \frac{2}{n} \boldsymbol{\varTheta} \sum\limits_{i=1}^{n} \boldsymbol{z}_{i}\boldsymbol{\phi}_{i}^{\top}; \end{array} $$

(34)

$$ \begin{array}{@{}rcl@{}} \nabla^{2} {\varLambda}(\boldsymbol{X}) &=& \left( \frac{2}{n} \sum\limits_{i=1}^{n}\boldsymbol{\phi}_{i} \boldsymbol{\phi}_{i}^{\top} \right) \otimes \boldsymbol{\varTheta}. \end{array} $$

Here ⊗ denotes the Kronecker product. Since ∇²Λ(X) is symmetric, we yield $\rho (\nabla ^{2} {\varLambda }(\boldsymbol {X})) = \|\nabla ^{2} {\varLambda }(\boldsymbol {X})\|_{2} = \|\left (\frac {2}{n} {\sum }_{i=1}^{n}\boldsymbol {\phi }_{i} \boldsymbol {\phi }_{i}^{\top } \right ) \otimes \boldsymbol {\varTheta }\|_{2} = \|\frac {2}{n} {\sum }_{i=1}^{n}\boldsymbol {\phi }_{i} \boldsymbol {\phi }_{i}^{\top }\|_{2} \|\boldsymbol {\varTheta }\|_{2}$. Thus, the proof is complete. □

Appendix B: Comparative Algorithms for Solving the Problem (2)

The Al-M method alternates between computing the variable X and Θ at each iteration. In particular, at iteration k, for fixed Θ, we need to compute X^k+ 1, an optimal solution to the following problem (see, e.g., [1, 44])

$$ \min \frac{1}{n} \sum\limits_{i=1}^{n} (\boldsymbol{z}_{i}-\boldsymbol{X} \boldsymbol{\phi}_{i})^{\top} \boldsymbol{\varTheta}^{k} (\boldsymbol{z}_{i}-\boldsymbol{X} \boldsymbol{\phi}_{i}) \text{ s.t. } \boldsymbol{X} \in \mathcal{X}. $$

(35)

Let us denote by Z (resp. Φ) a matrix in $\mathbb {R}^{m \times n}$ (resp. $\mathbb {R}^{d \times n}$) whose each column is a vector z_i (resp. ϕ_i); and define D^k := (ΦΦ^⊤)^(− 1/2)(ΦZ^⊤)(Θ^k)^(1/2). A reduced-rank regression estimate X^k+ 1 of (35) is given by (see [1])

$$ \boldsymbol{X}^{k+1} = \sum\limits_{t=1}^{r} {\varLambda}_{t} \left[(1/n) \boldsymbol{\Phi} \boldsymbol{\Phi}^{\top} \right]^{(-1/2)} \boldsymbol{u}_{t} \boldsymbol{v}_{t}^{\top} (\boldsymbol{\varTheta}^{k})^{(-1/2)}, $$

(36)

where the sequence {Λ_t} is the singular values of matrix D^k; {u_t} and {v_t} are the left-hand and right-hand singular vectors of D^k. For fixed X, the Al-M computes the point Θ^k+ 1 using (17) at X^k+ 1. Note that the Al-M method does not have any parameter.

The Al-GD method differs from the Al-M method by the fact that the Al-GD performs one iteration of gradient descent method for solving the problem (35) (see [9]). In particular, X^k+ 1 is computed as follows:

$$ \boldsymbol{X}^{k+1} = \text{Proj}_{\mathcal{X}}\left( \boldsymbol{X}^{k} + \frac{2\eta_{\boldsymbol{X}}}{n}\boldsymbol{\varTheta}^{k} \sum\limits_{i=1}^{n} (\boldsymbol{z}_{i}-\boldsymbol{X}^{k} \boldsymbol{\phi}_{i})\boldsymbol{\phi}_{i}^{\top} \right), $$

(35)

where the step size η_X is a tuning parameter. Al-GD computes the point Θ^k+ 1 using (17).

The J-GD method does not compute two variables alternately, but takes one gradient descent step in the joint variable (X,Θ) (see [9]). For estimating X^k+ 1, it is the same as (4), while the point Θ^k+ 1 is computed by using gradient descent method for (16) at the point (X^k, Θ^k) as follows:

$$ \boldsymbol{\varTheta}^{k+1} = \text{Proj}_{{\varOmega}}\left( \boldsymbol{\varTheta}^{k} + \eta_{\boldsymbol{\varTheta}} \boldsymbol{\Delta}^{k}\right), $$

(38)

where the step size η_Θ is a tuning parameter and

$$ \boldsymbol{\varDelta}^{k} = (\boldsymbol{\varTheta}^{k})^{(-1)} -\left[\frac{1}{n}\boldsymbol{\varTheta}^{k} \sum\limits_{i=1}^{n} (\boldsymbol{z}_{i}-\boldsymbol{X}^{k} \boldsymbol{\phi}_{i})(\boldsymbol{z}_{i}-\boldsymbol{X}^{k} \boldsymbol{\phi}_{i})^{\top} \right]. $$