## 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 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 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.

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