Prev: how to avoid NaN
Next: Cholesky decomposition
From: Anubhab on 29 Apr 2010 03:15
Hi,

I need to do a Cholesky decomposition of a matrix, Q. The matrix arises as a weight in computing energy in problems in hydrodynamic stability. Q is negative laplacian in 2d.
Q=-(D^2-k^2); (D is a Chebyshev differentiation matrix along y and k is a wavenumber along x obtained on Fourier transforming).
Energy= y'*Q*y.
I expect Q to be positive definite and turns out all its eigenvalues are >0. I have pasted the code below

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; clc
N=15; % number of discretization points

% Compute the Chebyshev differentiation matrix
x = cos(pi*(0:N)/N)';
c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
X = repmat(x,1,N+1);
dX = X-X';
D = (c*(1./c)')./(dX+(eye(N+1))); % off-diagonal entries
D = D - diag(sum(D,2)); % diagonal entries
D2 = D^2;

% Impose homogeneous boundary conditions
D2 = D2(2:N,2:N);
I = eye(size(D2)); % Identity matrix

k=1;
Q = -(D2-k^2*I);
all(eig(Q)>0)
chol(Q);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Now when N is small (~10) it works fine. As I increase N I get the following error

??? Error using ==> chol
Matrix must be positive definite.

whereas all the eigenvalues are still positive. Any suggestions why Cholesky decomposition isn't working?
Thanks

anubhab

p.s. The eigenvalues increase in magnitude with N. By N=20 the largest eigenvalues a little less than 10^4.
From: Joris de Wind on 29 Apr 2010 05:35
Did you check whether Q is symmetric??

If you read the help of chol (type help chol in the command window), you will notice that only the upper triangular part of Q is used to produce the Cholesky decomposition. Matlab assumes that the lower triangular part is the same as the upper triangular part.

If you look at your Q matrix, this is absolutely not the true. There are huge differences between Q and Q'. So probably you did make some typo somewhere in your code. The eigenvalues of Q may all be positive but Matlab is working on a totally different matrix (since it only uses the upper triangular part and assumes the lower triangular part is the same).

Just try to find where the construction of your Q matrix goes wrong.

Best, Joris





"Anubhab " <anubhab(a)jncasr.ac.in> wrote in message <hrbbnj$igh$1(a)fred.mathworks.com>...
> Hi,
>
> I need to do a Cholesky decomposition of a matrix, Q. The matrix arises as a weight in computing energy in problems in hydrodynamic stability. Q is negative laplacian in 2d.
> Q=-(D^2-k^2); (D is a Chebyshev differentiation matrix along y and k is a wavenumber along x obtained on Fourier transforming).
> Energy= y'*Q*y.
> I expect Q to be positive definite and turns out all its eigenvalues are >0. I have pasted the code below
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> clear all; clc
> N=15; % number of discretization points
>
> % Compute the Chebyshev differentiation matrix
> x = cos(pi*(0:N)/N)';
> c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
> X = repmat(x,1,N+1);
> dX = X-X';
> D = (c*(1./c)')./(dX+(eye(N+1))); % off-diagonal entries
> D = D - diag(sum(D,2)); % diagonal entries
> D2 = D^2;
>
> % Impose homogeneous boundary conditions
> D2 = D2(2:N,2:N);
> I = eye(size(D2)); % Identity matrix
>
> k=1;
> Q = -(D2-k^2*I);
> all(eig(Q)>0)
> chol(Q);
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> Now when N is small (~10) it works fine. As I increase N I get the following error
>
> ??? Error using ==> chol
> Matrix must be positive definite.
>
> whereas all the eigenvalues are still positive. Any suggestions why Cholesky decomposition isn't working?
> Thanks
>
> anubhab
>
> p.s. The eigenvalues increase in magnitude with N. By N=20 the largest eigenvalues a little less than 10^4.
From: John D'Errico on 29 Apr 2010 05:54
"Anubhab " <anubhab(a)jncasr.ac.in> wrote in message <hrbbnj$igh$1(a)fred.mathworks.com>...
> Hi,
>
> I need to do a Cholesky decomposition of a matrix, Q. The matrix arises as a weight in computing energy in problems in hydrodynamic stability. Q is negative laplacian in 2d.
> Q=-(D^2-k^2); (D is a Chebyshev differentiation matrix along y and k is a wavenumber along x obtained on Fourier transforming).
> Energy= y'*Q*y.
> I expect Q to be positive definite and turns out all its eigenvalues are >0. I have pasted the code below
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> clear all; clc
> N=15; % number of discretization points
>
> % Compute the Chebyshev differentiation matrix
> x = cos(pi*(0:N)/N)';
> c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
> X = repmat(x,1,N+1);
> dX = X-X';
> D = (c*(1./c)')./(dX+(eye(N+1))); % off-diagonal entries
> D = D - diag(sum(D,2)); % diagonal entries
> D2 = D^2;
>
> % Impose homogeneous boundary conditions
> D2 = D2(2:N,2:N);
> I = eye(size(D2)); % Identity matrix
>
> k=1;
> Q = -(D2-k^2*I);
> all(eig(Q)>0)
> chol(Q);
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> Now when N is small (~10) it works fine. As I increase N I get the following error
>
> ??? Error using ==> chol
> Matrix must be positive definite.
>
> whereas all the eigenvalues are still positive. Any suggestions why Cholesky decomposition isn't working?
> Thanks
>
> anubhab
>
> p.s. The eigenvalues increase in magnitude with N. By N=20 the largest eigenvalues a little less than 10^4.

A matrix that is THEORETICALLY positive definite can
still be numerically singular.

What matters here is not only the smallest eigenvalue,
but how large they get. Remember, you are working in
only about 16 digits of floating point arithmetic.

Finally, that last test, all(eig(Q)) > 0, is foolish.

(1+sqrt(-1)) > 0
ans =
1

Be very careful when you try to compare numbers
that may have an imaginary part. MATLAB compares
only the real parts!

John
From: Anubhab on 29 Apr 2010 09:58
Thanks for the comments.

@Joris

It is weird that I didn't notice what I was computing was not a symmetric matrix. I will go back and check where I possibly went wrong.

@John

Yeah, I guess it was probably not a very judicious thing to naively try the Cholesky decomposition when my eigenvalues were getting so large.
 | 
Pages: 1
Prev: how to avoid NaN
Next: Cholesky decomposition