3
$\begingroup$

I have 2 pixels with size 1x3 called $A$ and $B$ and I have to compute the following equation:

$$ A^T *(\Sigma+ I_3*\lambda)^{-1}*B $$

where $\Sigma$ is the covariance matrix (3x3) between vectors $A$ and $B$.

$I_3$ is the 3x3 identity matrix.

$\lambda$ is a constant (therefore, the matrix is not singular).

At the moment, i'm computing the inverse of the $\Sigma +I_3*\lambda$ using the Gauss-Jordan elimination.

I wanna know if there is a trick to compute this equation without computing the inverse. I'm also limited in memory so the Gauss-Jordan elimination is not a really good solution. I also tried to compute straight the inverse using the rule of Sarrus but the result was not enought accurate.

My aim is to resolve this equation with the highest speed and the minimum memory space.

EDIT:

Anyone knows a fast and good way to inverse a 3x3 symmetric matrix ?

EDIT 2:

I'm thinking about making a Cholesky decomposition of my matrix but after that, I don't understand how to compute the inverse of $(\Sigma +I_3*\lambda)$ from Cholesky matrix.

  • 0
    Hmm, set aside the question of whether Gauss-Jordan elimination is a good method, are you saying that your memory is so limited that you cannot even store a $3\times3$ matrix? Also, I am not an expert in numerical linear algebra, but computer algorithms usually have space-time tradeoffs. That you need something "with the highest speed and the minimum memory space" is perhaps asking for too much .... Anyway, if those matrix algorithms seem too complicated to you, I think you may try the "conjugate gradient method", which is fairly efficient in terms of speed, storage and ease of implementation.2012-11-15
  • 0
    @user1551 If i'm limited in memory is because I'm using this in a kernel function in graphic card. But I guess you're right about the trade off, that's why I'm trying different methods. About the "conjugate gradient method", it is usefull for equation Ax=b right ? How can I apply it for my equation ?2012-11-15
  • 0
    $x=(\Sigma+ I_3*\lambda)^{-1}*B \Leftrightarrow (\Sigma+ I_3*\lambda)x=B$.2012-11-15
  • 0
    @user1551 it works quite well, faster and less memory. However, the precondioned version is slower because it can already be solved in a few iterations. I will maybe try also using Cholesky to solve the same equation. I don't know which one would be the more efficient for my problem.2012-11-16
  • 0
    Hmm, if numerical stability is a concern, Cholesky decomposition seems to be inferior to QR factorization. See pp.17-18 of [these slides](http://www.ee.ucla.edu/~vandenbe/103/lectures/qr.pdf), for instance. If numerical stability is not your concern, actually you may just employ [Cramer's rule](http://en.wikipedia.org/wiki/Cramer%27s_rule).2012-11-16
  • 0
    Someone has posted on this forum [his MATLAB code for QR factorization](http://math.stackexchange.com/questions/238591/qr-decomposition-with-column-pivoting-errror/238785#238785) a few hours ago. You may go take a look.2012-11-16
  • 0
    @user1551 After using Cramer's rule, I have sometimes error like 10 times what I have with a cg. Is this can be du to numerical instability ? The results with Cramer's rule leads my algorithm to fail. I don't know if I'll have time to implement QR factorization but looks really interesting.2012-11-20
  • 0
    @user1551 Btw, you should put your comment as an answer like that I can accept it.2012-11-21
  • 0
    Thanks, but I don't think I have given any answer, so I'll pass. :-D2012-11-21

1 Answers 1

1

Since the matrix is symmetric, its inverse is also symmetricת we'll use that and solve.

Let this be the matrix:

[m11  m12  m13]   [m12  m22  m23]   [m13  m23  m33]   

Its determinant is:

D = m11 * (m33 * m22 - m23^2) - m21 * (m33 * m12 - m23 * m13) + m13 * (m23 * m12 - m22 * m13) 

Assuming it is non zero.

The members of the inverse:

a11 = m33 * m22 - m23^2   a12 = m13 * m23 - m33 * m12   a13 = m12 * m23 - m13 * m22    a22 = m33 * m11 - m13^2   a23 = m12 * m13 - m11 * m23  a33 = m11 * m22 - m12^2 

Divide them all by the determinant which is now given by:

D = (m11 * a11) + (m21 * a21) + (m13 * a13) 
  • 0
    m21 and a21 in the last D computation should presumably be m12 and a12 in the symmetrical notation.2018-03-12