I need to optimise the following scalar function with respect to a matrix $S$. $$ f(S) = \boldsymbol{y}^{T}\boldsymbol{X}w - \boldsymbol{1}_{n}^{T} \exp \left\{ \boldsymbol{X}w + \frac{1}{2} \operatorname{diag}(\boldsymbol{X}S\boldsymbol{X}^{T}) \right\} - \frac{1}{2}(w - \boldsymbol{\mu})^{T}\Sigma^{-1}(w-\boldsymbol{\mu}) - \frac{1}{2} \operatorname{tr}(\Sigma^{-1}S) + \frac{1}{2} \log \left| S \right| -\frac{1}{2}\log \left| \Sigma \right| $$

$\boldsymbol{y}: n \times1 \text{ column vector}$

$\boldsymbol{X}: n \times p \text{ matrix}$

$w, \, \boldsymbol{\mu}: p \times 1 \text{ column vectors}$

$S, \, \Sigma: p \times p \text{ semi-positive definite matrices}$

$\exp: \text{element-wise exponentiation}$

$\boldsymbol{1}_{n}: n \times 1 \text{ column vector with only ones}$

$\operatorname{diag}: \text{vectorizing the diagonal elements of the matrix inside}$

The Newton-Raphson method encompasses the optimisation of a scalar function with respect to a vector, resulting in computing the Jacobian matrix and then iteratively performing line search. However, here it's not a vector that is involved but rather a matrix. Is the Newton-Raphson method implemented equivalently, or is there something I should modify? (If there are other methods I can use, please tell me. It doesn't necessarily have to be Newton-Raphson.)


To save some typing, define
$$\eqalign{ E &= {\rm Diag}\Bigg(\exp\Big(Xw+\frac{1}{2}{\rm diag}(XSX^T)\Big)\Bigg) \cr M &= \Sigma^{-1} \cr }$$ Next, collect all of the pieces of the objective function which do not depend on $S$ into a constant term $K$, and write a simplified version $$\eqalign{ f &= K + \frac{1}{2}\Big({\rm tr}(\log(S)) - M:S\Big) - 1:E \cr }$$ Per your previous question, you already know the gradient of the final $E$-term. The gradient of the other 2 terms are simple. Thus the gradient of the overall function is $$\eqalign{ G &= \frac{\partial f}{\partial S} &= \frac{1}{2}\Big(S^{-1}-M-X^TEX\Big) \cr }$$ Any method based on the Hessian {Newton-Raphson, Quasi-Newton, Levenberg-Marquardt} is going to be impractical, since the Hessian will be a $(p\times p\times p\times p)\,$ 4th-order tensor --- or a $(p^2\times p^2)$ matrix if you apply vectorization tricks.

Also, finding a closed-form expression for the Hessian will be quite challenging. So that leaves you with gradient-based methods, which is okay.

Barzilai-Borwein is particularly easy to adapt to the situation where the working variables are matrices instead of vectors. The iteration is $$\eqalign{ dS_{k} &= S_{k} - S_{k-1} \cr dG_{k} &= G_{k} - G_{k-1} \cr \lambda_{k} &= \frac{dS_{k}:dG_{k}}{dG_{k}:dG_{k}} \cr S_{k+1} &= S_{k} - \lambda_{k} G_{k} \cr }$$ In place of the scalar product of 2 vectors, e.g. $(ds^Tdg),\,$ in the typical vector-based method, the above iteration substitutes the Frobenius product of the corresponding matrices, $(dS:dG)$.

The Barzilai-Borwein method works just fine for most problems, but may require an unconventional, non-monotone line search for "hard" problems. Google for "Raydan Barzilai" for details.

  • $\begingroup$ I'm actually optimising both $w$ and $S$. Would it be okay to combine the Newton-Raphson method using Hessian for $w$ and Barzilai-Borwein method for $S$? $\endgroup$ – Daeyoung Lim Dec 14 '15 at 7:02
  • $\begingroup$ And is it correct that I have to initialise 2 of $S$ ($S_{0}, S_{1}$) in order to calculate the first $dS_{1}$? $\endgroup$ – Daeyoung Lim Dec 14 '15 at 9:30
  • $\begingroup$ How can I adjust for large values resulting from the exponentiation in the gradient of $S$, $\frac{\partial f}{\partial S}$? I've thought of taking the logarithm but it's not appropriate here since it has to be multiplied by $X^{T}$ and $X$... Is there any way I can circumvent this computationally? $\endgroup$ – Daeyoung Lim Dec 14 '15 at 12:29
  • 1
    $\begingroup$ Unless the function is separable, you cannot use Newton-Raphson on $w\,$ and independently use Barzilai-Borwein on $S$. Unfortunately, your function is not separable, since both variables appear in the exponential term. You'll need to reformulate your function in terms of a single independent variable $V$ (which somehow incorporates both $w$ and $S$), and find the gradient with respect to this new variable. $\endgroup$ – lynn Dec 14 '15 at 17:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.