I am trying to implement an algorithm from a paper which makes use a QR factorization of a real matrix $A$ as a means of one of forming the Cholesky factor of $A^T A$ without explicitly forming $A^T A$ (which is generally ill-conditioned).
In particular:
We have a covariance matrix $\mathbf{P}_x = \sum_{i=0}^{2L} w_i \mathbf{x}_i \mathbf{x}_i^T$ and we want to find the upper Cholesky factorization $\mathbf{S}_x^T \mathbf{S}_x = \mathbf{P}_x$. The first weight, $w_0$, may be negative but the rest are positive.
The technique advocated in the paper is a QR decomposition followed by a Cholesky update (or downdate) for the negative weight. That is, let $\mathbf{P}_x = \mathbf{P}_0 + \mathbf{P}_i$ where $\mathbf{P}_0 = w_0 \mathbf{x}_0\mathbf{x}_0^T$ and $\mathbf{P}_i = \sum_{i=1}^{2L} w_i \mathbf{x}_i \mathbf{x}_i^T$. So if $\mathbf{S_0}^T \mathbf{S_0} = \mathbf{P}_0$, then (using MATLAB's syntax) $\mathbf{S}_x = \operatorname{cholupdate}(\mathbf{S}_0, \sqrt{\|w_0\|}\mathbf{x}_0, \operatorname{sign}(w_0))$.
The trick comes from forming $\mathbf{S}_0$ without the Cholesky decomposition. Define $\mathbf{A} = \sum_{i=1}^{2L} \sqrt{w_i} \mathbf{x}_i$, and let $\mathbf{A} = \mathbf{Q}\mathbf{R}$. It can be seen that $\mathbf{A}^T \mathbf{A} = \mathbf{S}_0^T \mathbf{S}_0 = (\mathbf{Q} \mathbf{R})^T \mathbf{Q} \mathbf{R} = \mathbf{R}^T \mathbf{Q}^T \mathbf{Q} \mathbf{R} = \mathbf{R}^T \mathbf{R}$, and so I would expect $\mathbf{R} = \mathbf{S}_0$.
However, my MATLAB implementation shows that $\mathbf{R} \neq \mathbf{S}_0$. Am I incorrect in assuming that the Cholesky factorization is unique, or is there a subtlety that I am missing?
Here is an example MATLAB script which demonstrates the issue:
% Weights
w = [-3, 1, 1, 1, 1, 1, 1];
% Generate a whole bunch of random vectors
A = rand(3, 7)
% Calculate the full P Matrix
Px = zeros(3,3);
for i = 1:7
Px = Px + w(i) * A(:,i) * A(:,i).';
end
% Calculate part of QR
A_dash = A(:, 2:7);
B = triu(qr(sqrt(w(2)) * A_dash.', 0));
% Cholupdate it
S = cholupdate( B(1:3, 1:3), sqrt(abs(w(1))) * A(:,1), '-')
P = S.' * S
% DisplayCholesky of full covariance
Px
Sx = chol(Px)
The output from a particular run was:
A =
0.1788 0.5985 0.6999 0.0688 0.6544 0.7184 0.3251
0.4229 0.4709 0.6385 0.3196 0.4076 0.9686 0.1056
0.0942 0.6959 0.0336 0.5309 0.8200 0.5313 0.6110
S =
-1.3442 -1.1314 -1.1479
0 -0.1750 -0.0802
0 0 -0.8629
P =
1.8070 1.5209 1.5430
1.5209 1.3107 1.3128
1.5430 1.3128 2.0686
Px =
1.8070 1.5209 1.5430
1.5209 1.3107 1.3128
1.5430 1.3128 2.0686
Sx =
1.3442 1.1314 1.1479
0 0.1750 0.0802
0 0 0.8629
So indeed, the covariance matrices are the same, but the Cholesky factors differ by a sign. I have a couple of questions:
- Is this the result of a QR decomposition not being uniquely defined? I thought the Cholesky Decomposition was uniquely defined!
- If so, what routine or method (MATLAB, C/C++, LAPACK) will obtain the same Cholesky factor?
- Or, is there a better algorithm altogether?