sqrtpsdm

Produces a "square-root" of a positive semi-definite matrix, S, such that:

 S == sqS * sqS.'

When S has n eigenvalues of 0, the square-root will not be unique, and n of the columns of sqS will be zeros. (Usually, the first n colums will contain the zeros, but this is not true in certain pathological cases). The function also returns the number of zeros it finds.

sqS = sqrtpsdm(S)
[sqS, nz] = sqrtpsdm(S)

Note that this will still be a matrix square-root: sqS * sqS.' == S.

A different meaning of "matrix square root" would be C in C * C == S; pleae note that this is not the form intended here.

Inputs

S	A positive, semi-definite matrix
type	Use 'L' for a lower-triangular square root and 'U' for an upper-triangular square root. When using 'L', and when the matrix is positive-definite, `sqS` will be the lower Cholesky factor (and, therefore, the transpose of the upper Cholesky factor).

A positive, semi-definite matrix

type

Use 'L' for a lower-triangular square root and 'U' for an upper-triangular square root.

When using 'L', and when the matrix is positive-definite, sqS will be the lower Cholesky factor (and, therefore, the transpose of the upper Cholesky factor).

Outputs

sqS	A matrix square root of `S`
n	Number of 0 eigenvalues of `S`

Example: Nominal Case

Make a positive definite matrix, find the sqrt, and reconstruct.

P = randcov(4)
sqP = sqrtpsdm(P);
sqP * sqP.'

P =
    0.4875   -0.1141   -0.3501   -0.0852
   -0.1141    0.4080    0.0600    0.3106
   -0.3501    0.0600    0.6409   -0.0734
   -0.0852    0.3106   -0.0734    0.5711
ans =
    0.4875   -0.1141   -0.3501   -0.0852
   -0.1141    0.4080    0.0600    0.3106
   -0.3501    0.0600    0.6409   -0.0734
   -0.0852    0.3106   -0.0734    0.5711

Example: A 0 Eigenvalue

Make a positive, semi-definite matrix, find the sqrt, and reconstruct. Also, have it report the number of zeros found.

P = diag([1 2 3 0])
[sqP, nz] = sqrtpsdm(P);
sqP * sqP.'
nz

P =
     1     0     0     0
     0     2     0     0
     0     0     3     0
     0     0     0     0
ans =
    1.0000         0         0         0
         0    2.0000         0         0
         0         0    3.0000         0
         0         0         0         0
nz =
     1

Show that the Cholesky decomposition won't work because we'd need a fully positive-definite matrix.

try
    chol(P)
catch err
    fprintf('Caught an error from chol: %s\n', err.message);
end

Caught an error from chol: Matrix must be positive definite.

Example: Another Bad Eigenvalue

That was an easy case. Let's try rotating the matrix so the problematic eigenvalue isn't obvious.

theta = 2*pi*rand();
R     = [ cos(theta), sin(theta); ...
         -sin(theta), cos(theta)];
A = R * [-eps 0; 0 1] * R.' % This has 1 "bad" eigenvalue.
[sqA, nz] = sqrtpsdm(A);
A * A.'
nz

A =
    0.8494   -0.3576
   -0.3576    0.1506
ans =
    0.8494   -0.3576
   -0.3576    0.1506
nz =
     1

Example: Other Forms

By default, sqrtpsdm uses a UDU factorization. We can instead use LDL if we desire.

P = randcov(3);
sqP_upper = sqrtpsdm(P)
sqP_lower = sqrtpsdm(P, 'lower')

sqP_upper =
    0.7454    0.0254    0.1062
         0    0.8931    0.1429
         0         0    0.8930
sqP_lower =
    0.7533         0         0
    0.0502    0.9030         0
    0.1259    0.1343    0.8738

They're still all the same.

P
sqP_upper * sqP_upper.'
sqP_lower * sqP_lower.'

P =
    0.5675    0.0378    0.0949
    0.0378    0.8180    0.1276
    0.0949    0.1276    0.7974
ans =
    0.5675    0.0378    0.0949
    0.0378    0.8180    0.1276
    0.0949    0.1276    0.7974
ans =
    0.5675    0.0378    0.0949
    0.0378    0.8180    0.1276
    0.0949    0.1276    0.7974

Docs