2
$\begingroup$

I'm exploring different techniques to efficiently solve some matrix equations. My situation is that I have a matrix $\textbf{H} = \textbf{J}^T \textbf{J}$, where $\textbf{J}$ is a matrix with no special qualities. From this, I first need to calculate $\mathrm{tr}(\textbf{H}^{-1})$, and then I solve $(\textbf{H}+\lambda \textbf{I})\vec{\delta} = \vec{g}$, where $\vec{g}$ is known, repeatedly using several different values of $\lambda$. (If it helps, all values are real here).

If I understand correctly, $\textbf{H}$ is guaranteed symmetric and positive-semidefinite. Therefore, I know that I can use a cholesky or LDLT decomposition for each of the calculations, however I'd rather not perform this many decompositions. Am I right in saying that, because $\textbf{H}$ is symmetric and real, that I can perform an eigendecomposition on it to get $\textbf{H} = \textbf{Q}\mathbf{\Lambda}\textbf{Q}^{-1}$ for all $\textbf{H}$? And if so, is the following correct?

$$\mathrm{tr}(\textbf{H}^{-1}) = \sum_i \frac{1}{\mathbf{\Lambda}_i}$$ $$\textbf{H}+\lambda \textbf{I} = \textbf{Q}(\mathbf{\Lambda}+\lambda\textbf{I})\textbf{Q}^{-1}$$ $$\vec{\delta} = \vec{g}\textbf{Q}(\mathbf{\Lambda}+\lambda\textbf{I})^{-1}\textbf{Q}^{-1}$$

Lastly, there is a variation of this algorithm using $\textrm{diag}(\mathbf{H})$ instead of $\mathbf{I}$. Can I perform a similar eigendecomposition on $\textbf{H}+\lambda \textrm{diag}(\textbf{H})$ to get a constant factor $\textbf{Q}$ and the eigenvalues for all values of $\lambda$ with a single decomposition?

Whilst I realise that eigendecompositions are much slower than an LDLT decomposition, I believe that the above transformation would allow me to perform significantly less decompositions, and thus hopefully improve the efficiency of this algorithm.

Improve this question
$\endgroup$
2
  • $\begingroup$ Without claiming to have been careful, this looks OK to me, except I don't think the part beginning "Lastly" works. $\endgroup$
    – Steve Huntsman
    Jul 23, 2013 at 13:34
  • $\begingroup$ Ok great :) Yeah I didn't think so, was just wondering if there was some property I had overlooked that could be applied there. $\endgroup$
    – sourtin
    Jul 23, 2013 at 13:43

2 Answers 2

Reset to default
1
$\begingroup$

If $D = diag(H)$ is positive, let $K = D^{-1/2} H D^{-1/2}$ which again is positive semidefinite. Then $(H + \lambda D)^{-1} = D^{-1/2} (K + \lambda I)^{-1} D^{-1/2}$, so if you do the eigendecomposition $K = Q \Lambda Q'$ with $Q$ orthogonal you get $(H+\lambda D)^{-1} = (D^{-1/2} Q) (\Lambda + \lambda I)^{-1} (D^{-1/2} Q)'$.

Improve this answer
$\endgroup$
1
  • $\begingroup$ That's perfect, thanks! I deleted my last comment because I just realised how easy $\textrm{tr}(\mathbf{H}^{-1})$ is to find using $\mathbf{K}$. $\endgroup$
    – sourtin
    Jul 23, 2013 at 17:41
1
$\begingroup$

I would advise against computing $H$. Rather, compute SVD of $J$, i.e., $J = U \Sigma V^T$, so $$H = J^T J = V \Sigma^T U^T U \Sigma V^T = V \Sigma^T \Sigma V^T.$$ If $J$ is square or tall, then $\Lambda = \Sigma^2$. Otherwise, you get additional zeroes (which do not affect the trace, but in this case $H$ would be singular).

Improve this answer
$\endgroup$
5
  • $\begingroup$ That's a really interesting approach, but I think it would be less efficient than using an eigendecomp on $H$. First, as $J$ may have more rows than $H$ (Currently I build $H$ up iteratively instead of storing the whole of $J$) it would need more memory. Second, my matrix library says it's slower for SVD than ED. I will try it though, thanks for the answer! $\endgroup$
    – sourtin
    Jul 23, 2013 at 17:25
  • 1
    $\begingroup$ First, you can use reduced SVD, or some other variant, to save on computation. Note that you don't need $U$ at all and, if it's just the trace you're after, you need to compute only $\Sigma$, which is of the same order as $H$. This means no additional use of memory (just keep changing $J$ until you get a diagonal). Second, computing $H$ takes time and introduces errors. One of the main applications of SVD is finding the eigenvalues of $H = J^TJ$ without actually computing $H$. For example, see the end of the first answer here. $\endgroup$
    – Vedran Šego
    Jul 23, 2013 at 17:36
  • $\begingroup$ Looking into using SVD further, I can combine both yours and Robert's method in order to use the $\textrm{diag}(\mathbf{H})$ variation by performing SVD on $\mathbf{JD}^{-1/2}$ instead of on $\mathbf{J}$. I plan on testing all these different methods to see which performs better, and the improved accuracy may reduce the number of iterations and so make up for being slower. It'll be interesting to see the difference, as $H=J^TJ$ is itself a (good) approximation, but perhaps the better handling of cases such as the Läuchli matrix will be significant. $\endgroup$
    – sourtin
    Jul 23, 2013 at 19:11
  • $\begingroup$ Sorry that I haven't accepted your answer, it has been very valuable, but Robert more directly addressed my question. If yours performs better on my benchmarks though then I will come back and accept yours. $\endgroup$
    – sourtin
    Jul 23, 2013 at 19:15
  • $\begingroup$ Don't worry about that. I'm glad I've been able to help and Robert did give a good answer. $\endgroup$
    – Vedran Šego
    Jul 23, 2013 at 19:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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