How does Cholupdate work?

I have made some MATLAB code for Cholupdate algorithm. The problem is that some times it works, some times not.

function [L] = cholupdate2(L, x, OP)
  cholupdate(L, x, OP) % OP can be "+" = rank_one_update or "-" = rank_one_downdate
  F = 1;
  if(strcmp(OP, '-') == 1)
    F = -1; % Rank_one_downdate
  end
  n = length(x);
  for k = 1:n
      a = L(k,k);
      b = x(k);
      r = sqrt(a*a + F * b*b);
      c = r / a;
      s = b / a;
      L(k, k) = -r;
      if k < n
          for i = k+1:n
            L(k, i) = -(L(k,i) + F * s * x(i)) / c;
            x(i) = -c * x(i) - s * L(k, i);
          end
      end
  end
end

If I try an example

A = [-3.2462,  2.1783,      2.4925,           0;
           0, -3.3096,     -1.4982,           0;
           0,       0,     -5.8431,     -1.0139;
           0,       0,           0,     -4.8889];

x =  [-3.0105,
       2.2319,
       2.5538,
           0];

cholupdate2(A, x, '+')  %// or '-' (see the MATLAB code in cholupdate function above)

Then I get the output.

ans =

  -4.4273   3.1148   3.5641        0
        0  -3.3132  -1.5049        0
        0        0  -5.8441  -1.0137
        0        0        0  -4.8889

ans =

  -4.4273   3.1148   3.5641        0
        0  -3.3132  -1.5049        0
        0        0  -5.8441  -1.0137
        0        0        0  -4.8889

The first matrix is from MATLAB/GNU Octave, the second matrix is from my cholupdate2 function.

But if I change -3.3096 to 3.3096 inside the A matrix.

Now I'm getting the output. As you can see, they are not 100% exactly the same.

ans =

  -4.4273   3.1148   3.5641        0
        0   3.3132  -1.4882        0
        0        0  -5.8483  -1.0130
        0        0        0  -4.8891

ans =

  -4.4273   3.1148   3.5641        0
        0  -3.3132   1.4882        0
        0        0  -5.8483  -1.0130
        0        0        0  -4.8891

I took the source code from Wikipedia and modify so it will return not-transpose back, which is do inside the Wikipedia code. I also added the rank-one downdate functionality.

So my question for you is:

If I got a matrix A and I want to update it with cholesky update. I need to add

$$\tilde A = A + x*x^T$$

Is that correct? What can $x$ be then?

Solution 1:

Consider the SPD matrix $\mathbf{A} = \mathbf{U}^T \mathbf{U}$. We are interested in the Cholesky decomposition of its rank-one update: $\mathbf{A}_1= \mathbf{A}+\mathbf{x}\mathbf{x}^T = \mathbf{U}_1^T \mathbf{U}_1$.

The cholupdate function computes $\mathbf{U}_1$ from the knowledge of $\mathbf{U}$ and $\mathbf{x}$

The idea is to consider the extended $(n+1)\times n$ matrix $\mathbf{B}= \begin{bmatrix} \mathbf{x}^T \\ \mathbf{U} \end{bmatrix} $ and find the (orthogonal) Givens rotations such that $ \mathbf{C} = \mathbf{G}_n \mathbf{G}_1 \mathbf{B}= \begin{bmatrix} \mathbf{0}^T \\ \mathbf{U}_1 \end{bmatrix} $ where $\mathbf{U}_1$ is upper triangular.

Because the Givens rotations are orthogonal, we see that $ \mathbf{C}^T\mathbf{C}= \mathbf{B}^T\mathbf{B}= \mathbf{U}^T\mathbf{U}+\mathbf{x}\mathbf{x}^T = \mathbf{U}_1^T\mathbf{U}_1$

Solution 2:

Your decomposition is as valid as the one from Matlab/GNU Octave because there is some degree of freedom in the row sign of $\mathbf{U}$. On your example, the Given matrix transforms

B=
   -3.0105    2.2319    2.5538         0
   -3.2462    2.1783    2.4925         0
         0   -3.3096   -1.4982         0
         0         0   -5.8431   -1.0139
         0         0         0   -4.8889

into

G1*B=
         0    0.1553    0.1776         0
   -4.4273    3.1148    3.5641         0
         0   -3.3096   -1.4982         0
         0         0   -5.8431   -1.0139
         0         0         0   -4.8889

then

G2*G1*B =
         0         0    0.1072         0
   -4.4273    3.1148    3.5641         0
         0   -3.3132   -1.5049         0
         0         0   -5.8431   -1.0139
         0         0         0   -4.8889

and finally

G3*G2*G1*B =

         0         0   -0.0000   -0.0186
   -4.4273    3.1148    3.5641         0
         0   -3.3132   -1.5049         0
         0         0   -5.8441   -1.0137
         0         0         0   -4.8889