Take the 2-minute tour ×
MathOverflow is a question and answer site for professional mathematicians. It's 100% free, no registration required.

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.

share|improve this question
    
Without claiming to have been careful, this looks OK to me, except I don't think the part beginning "Lastly" works. –  Steve Huntsman Jul 23 '13 at 13:34
    
Ok great :) Yeah I didn't think so, was just wondering if there was some property I had overlooked that could be applied there. –  nitrous Jul 23 '13 at 13:43
add comment

2 Answers

up vote 1 down vote accepted

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)'$.

share|improve this answer
    
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}$. –  nitrous Jul 23 '13 at 17:41
add comment

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).

share|improve this answer
    
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! –  nitrous Jul 23 '13 at 17:25
1  
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. –  Vedran Šego Jul 23 '13 at 17:36
    
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. –  nitrous Jul 23 '13 at 19:11
    
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. –  nitrous Jul 23 '13 at 19:15
    
Don't worry about that. I'm glad I've been able to help and Robert did give a good answer. –  Vedran Šego Jul 23 '13 at 19:40
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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