Eig function and chol function returning wrong matrix factorization

Hello, I simply tried to decompose a hermitian matrix C_z as [V_z,D_z] = eig(C_z) and compute L_z = V_z*sqrt(D_z).

I could also have done L_z = chol(C_z).

Surprisingly, the result is far from accurate. I got very large errors using eig and chol : L_z*L_z' is far from C_z.

In comparison, svd function is performing quite well.

With SVD: ||C_z - L_z.L_z^H|| = 1.796903e-13

With EIG: ||C_z - L_z.L_z^H|| = 1.579992e+02

With CHOL: ||C_z - L_z.L_z^H|| = 2.015977e+02

C_z is well conditioned, I tried to figure out where this came from but it seems that there must be a bug in the two functions (eig and chol). How do you explain this error ?

The code is here :

%--------------------------------------------------------------------------
% Compute pos
%--------------------------------------------------------------------------
Mx = 4;
My = 4;
x = (1:Mx)-1;
y = (1:My)-1;
pos = reshape([repmat(x,My,1),repmat(reshape(y,[],1),1,Mx)],My*Mx,2);
%--------------------------------------------------------------------------
% Compute arr
%--------------------------------------------------------------------------
x = -sqrt(2)/2:sqrt(2)/2;
y = -sqrt(2)/2:sqrt(2)/2;
Qx = length(x);
Qy = length(y);
arr = reshape([repmat(x,Qy,1),repmat(reshape(y,[],1),1,Qx)],Qy*Qx,2);
%--------------------------------------------------------------------------
% Compute x
%--------------------------------------------------------------------------
x = zeros(Qx,Qy);
x(round(Qy/4),(round(Qx/4+1))) = 10^(-10/10);
x(round(Qy/4),(round(Qx*3/4))) = 10^(-5/10);
x(round(Qy*3/4),(round(Qx/4+1))) = 10^(0/10);
x(round(Qy*3/4),(round(Qx*3/4))) = 10^(5/10);
x = x(:);
%--------------------------------------------------------------------------
% Compute A
%--------------------------------------------------------------------------
A = exp(pos*arr.');
%--------------------------------------------------------------------------
% Compute C_n
%--------------------------------------------------------------------------
C_n = 1e2 * eye(Mx*My);
%--------------------------------------------------------------------------
% Compute C_z and compare with L_z.L_z^H
%--------------------------------------------------------------------------
C_z = A*diag(x)*A'+C_n;
%%%%%
[V_z,D_z] = eig(C_z);
L_z_ = V_z*sqrt(D_z);
%%%%%%
[U_z,S_z,V_z] = svd(C_z);
L_z = U_z*sqrt(S_z);
%%%%%%
R = chol(C_z);
%%%%%%
fprintf(['With SVD: ||C_z - L_z.L_z^H|| = %i\n' ...
    'With EIG: ||C_z - L_z.L_z^H|| = %i\n' ...
    'With CHOL: ||C_z - L_z.L_z^H|| = %i\n\n'], norm(C_z - L_z*L_z'), norm(C_z - L_z_*L_z_'), norm(C_z - R*R'));

With SVD: ||C_z - L_z.L_z^H|| = 1.935988e-13
With EIG: ||C_z - L_z.L_z^H|| = 1.318081e+02
With CHOL: ||C_z - L_z.L_z^H|| = 2.015977e+02

Answer 1

[Bruno Luong]

You make two mistakes, please see comments (where "Bruno" appears) in this corrected code

%--------------------------------------------------------------------------
% Compute pos
%--------------------------------------------------------------------------
Mx = 4;
My = 4;
x = (1:Mx)-1;
y = (1:My)-1;
pos = reshape([repmat(x,My,1),repmat(reshape(y,[],1),1,Mx)],My*Mx,2);
%--------------------------------------------------------------------------
% Compute arr
%--------------------------------------------------------------------------
x = -sqrt(2)/2:sqrt(2)/2;
y = -sqrt(2)/2:sqrt(2)/2;
Qx = length(x);
Qy = length(y);
arr = reshape([repmat(x,Qy,1),repmat(reshape(y,[],1),1,Qx)],Qy*Qx,2);
%--------------------------------------------------------------------------
% Compute x
%--------------------------------------------------------------------------
x = zeros(Qx,Qy);
x(round(Qy/4),(round(Qx/4+1))) = 10^(-10/10);
x(round(Qy/4),(round(Qx*3/4))) = 10^(-5/10);
x(round(Qy*3/4),(round(Qx/4+1))) = 10^(0/10);
x(round(Qy*3/4),(round(Qx*3/4))) = 10^(5/10);
x = x(:);
%--------------------------------------------------------------------------
% Compute A
%--------------------------------------------------------------------------
A = exp(pos*arr.');
%--------------------------------------------------------------------------
% Compute C_n
%--------------------------------------------------------------------------
C_n = 1e2 * eye(Mx*My);
%--------------------------------------------------------------------------
% Compute C_z and compare with L_z.L_z^H
%--------------------------------------------------------------------------
C_z = A*diag(x)*A'+C_n;

% Change by Bruno, make C_z numerically Hermitian, otherwise eig won't
% return orthonormal V_z
C_z = 0.5*(C_z+C_z');

%%%%%
[V_z,D_z] = eig(C_z);
L_z_ = V_z*sqrt(D_z);
%%%%%%
[U_z,S_z,V_z] = svd(C_z);
L_z = U_z*sqrt(S_z);
%%%%%%
R = chol(C_z);
Lzchol = R'; % Bruno's comment: You make a mistake here
%%%%%%
fprintf(['With SVD: ||C_z - L_z.L_z^H|| = %i\n' ...
    'With EIG: ||C_z - L_z.L_z^H|| = %i\n' ...
    'With CHOL: ||C_z - L_z.L_z^H|| = %i\n\n'], norm(C_z - L_z*L_z'), norm(C_z - L_z_*L_z_'), norm(C_z - Lzchol*Lzchol'));

With SVD: ||C_z - L_z.L_z^H|| = 1.726597e-13
With EIG: ||C_z - L_z.L_z^H|| = 7.510780e-13
With CHOL: ||C_z - L_z.L_z^H|| = 3.181917e-14

[Cyril]

Thank you for the fast reply, but my concern was about the importance of the difference between svd end eig function results.

I think that, if such small numerical errors can lead to this big deviation, the problem must be adressed in the code of eig directly.

Here, all the code was simplied and the different values where changed compared to when a really faced the problem, and it persisted.

In addition, the matrix is (quite well) numericaly hermitian : norm(C_z - C_z') ~ 10^-15.

I suggest that if

C_z = 0.5*(C_z+C_z');

is necessary, it should be added in eig function.

[Bruno Luong]

"is necessary, it should be added in eig function."

NO not MATLAB fault, it's your mistake let me explain

Because the EIG factorization is this

C_z = V_z * D_z * inv(V_z)

This factorization is accurate. But you do not check that.

However what YOU assume is this

L_z*L_z' = V_z * D_z * V_z', but it is NOT C_z, compare with the above,

>> [V_z,D_z] = eig(C_z); % from C_z NOT symmetrized numerically
>> norm(V_z*D_z*inv(V_z)-C_z) % OK

ans =

   2.9983e-13

>> norm(V_z*D_z*V_z'-C_z) % WRONG

ans =

  157.9992

unless inv(V_z) = V_z'. (in other words V_z*V_z' = eye(n), we have othonormal eigen vectors)

If C_z is NOT numerical hermitian; EIG will not return orthonormal eigen vectors (they are not unique).

>> norm(eye(16)-V_z*V_z') % from C_z NOT symmetrized numerically

ans =

    1.5800

>> [V_z,D_z] = eig(0.5*(C_z+C_z'));
>> norm(eye(16)-V_z*V_z')

ans =

   1.3504e-15

Therefore your expectation of L_z_*L_z_' giving back C_z is wrong, unless you make C_z numerically Hermitian as I showed.