how do i solve the matrix equation
Show older comments
(X')*A*X=I where I is the unit matrix and A is positive definite and Hermitian
Answers (2)
I know this question is more than a decade old. For anyone trying to implement the top solution in MATLAB, note that there is a slight syntax typo there. The proposed solution I\sqrtm(A) evaluates to just sqrtm(A), but I believe it was intended to be sqrtm(A)\I to achieve the inverse 
. Else, the literal proposed solution only works in the trivial case where 
, as shown in Case 1.
. Else, the literal proposed solution only works in the trivial case where , as shown in Case 1.With that syntax correction in Case 2, the math holds true under the conjugate transpose 
. However, it is worth noting that this is not a unique solution. Because matrix multiplication is non-commutative, there are infinitely many valid matrices 
that satisfy this equation, as demonstrated in Case 3.
. However, it is worth noting that this is not a unique solution. Because matrix multiplication is non-commutative, there are infinitely many valid matrices that satisfy this equation, as demonstrated in Case 3.While Cholesky is the most common method taught in undergrad linear algebra class (for algorithmic efficiency), I believe there are several unexplored alternative methods to parameterize and generate the infinite solutions, especially if the matrix 
is close to being singular (See Case 4).
is close to being singular (See Case 4).Case 1: when is the identity matrix
is the identity matrix% Hermitian positive-definite matrix
A = eye(3)
tof = ishermitian(A) % returns logical 1 if matrix A is legitimately Hermitian
% Identity matrix (a.k.a unit matrix)
I = eye(size(A))
% Solution
X = I\sqrtm(A)
% Test: returns 0 if the result is correct
(X')*A*X - I
Case 2: is any arbitrary Hermitian positive-definite matrix (non-Identity) with 
is any arbitrary Hermitian positive-definite matrix (non-Identity) with % Hermitian positive-definite matrix
A = [2.5, 2.5, 0.5
2.5, 5.0, 1.5
0.5, 1.5, 2.0];
eig(A)
tof = ishermitian(A) % returns logical 1 if matrix A is legitimately Hermitian
% Identity matrix (a.k.a unit matrix)
I = eye(size(A))
% Solution
X = sqrtm(A)\I % or inv(sqrtm(A))
% Test: returns 0 (or near double-precision machine epsilon) if the result is numerically accurate
(X')*A*X - I
Case 3: using the inverse Cholesky factorization trick
% Cholesky decomposition (A = L*L')
L = chol(A, 'lower'); % same matrix A from Case 2
X_baseline = inv(L');
% generate a random unitary orthogonal matrix Q via QR decomposition
[Q, ~] = qr(rand(3)); % Q'*Q = I
% general solution
X = X_baseline*Q
% Test: returns 0 (or near double-precision machine epsilon) if the result is numerically accurate
(X')*A*X - I
Case 4: when is close to being singular
is close to being singularformat longg
lambda_max = 5 - 1e-10; % Largest eigenvalue
lambda_mid = 2 - 1e-11; % Middle eigenvalue
lambda_min = 1e-9; % Smallest eigenvalue
% diagonal eigenvalue matrix (Sigma)
Sigma = diag([lambda_max, lambda_mid, lambda_min]);
% construct a real Hermitian positive-definite matrix using QR decomposition trick
[Q, ~] = qr(rand(3));
A = Q*Sigma*Q'
cond(A)
eig(A)
tof = ishermitian(A) % returns logical 1 if matrix A is legitimately Hermitian
% Cholesky decomposition (A = L*L')
L = chol(A, 'lower'); % same matrix A from Case 2
X_baseline = inv(L');
% generate a random unitary orthogonal matrix Q via QR decomposition
[Q, ~] = qr(rand(3)); % Q'*Q = I
% general solution
X = X_baseline*Q
% Test: returns 0 (or near double-precision machine epsilon) if the result is numerically accurate
(X')*A*X - I
X=sqrtm(A)\I;
Categories
Find more on Linear Algebra in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
