Seemingly non-unique Cholesky factor via QR rectangularisation

Question

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?

Answer 1

The issue is that standard implementations of the QR decomposition do not require that $R$ has a positive diagonal. I recommend reading LAPACK Working Note 203, which introduces an alternate algorithm which guarantees a positive diagonal, namely xLARFP (e.g., zlarfp).

If you use a QR decomposition which results in an $R$ with a positive diagonal, your Cholesky factors will agree (up to floating point error). Traditionally, Cholesky factors are defined to have positive diagonals (as a result of taking square roots of positive entries), but obviously if $A=L L^H$, then $A=(-L)(-L)^H$ also holds. More generally, $A=(LD)(LD)^H$ for arbitrary diagonal $D$ whose diagonal entries lie on the complex unit circle.

Answer 2

Cholesky factorizations are only unique if the matrix $\mathbf{A}$ is Hermitian positive definite. Cholesky factorizations exist in the Hermitian positive semidefinite case, but then these factorizations are nonunique. Cholesky decompositions do not exist in all other cases.

The QR factorization of an full rank matrix is unique, if we restrict the diagonal of $\mathbf{R}$ to be positive.

Answer 3

From Wikipedia:

The Cholesky decomposition is unique: given a Hermitian, positive-definite matrix $A$, there is only one lower triangular matrix $L$ with strictly positive diagonal entries such that $A = LL^*$.

One of your "Cholesky factors" has a negative diagonal.