0

Færdig

# chol (~call) again

Pavel, for you in no case should underline the important function chol; unfortunately there are still some problems...

Please run the following few lines on the last version.

% *** ASSERTIONS *** (This has been incorporated into the documentation in R14SP3):

%

% CHOL expects its input matrix to be symmetric and only looks at the

% upper triangular portion of the matrix.

% CHOL function will return an error if it is only provided with a single

% output argument, and is also given a matrix that is not positive

% definite.

% CHOL function provides an optional second output argument "p" which

% is zero if the matrix is found to be positive definite. If the input

% matrix is not positive definite, then "p" will be a positive integer.

% *** ASSUMPTION ***

% Your implementation uses ?potr routine ...

% >>>>>>>>>> Let's test ...

%% A is scalar

% clc

mp.FollowMatlabNumericFormat(1) % just for convenience

A = mp(9) %#ok<*NOPTS>

R = chol(A) %#ok<*NASGU>

L = chol(A,'lower')

R = chol(A,'upper')

try [L,p] = chol(A,'lower'), catch, fprintf('mp [L,p] = chol(A,''lower'') BOOOM !!!\n\n'),...

[L,p] = chol(double(A),'lower'), fprintf('matlab ok\n\n'), end %#ok<*ASGLU>

try [P,p] = chol(A,'upper'), catch, fprintf('mp [L,p] = chol(A,''upper'') BOOOM !!!\n\n'), ...

[L,p] = chol(double(A),'upper'), fprintf('matlab ok\n\n'), end

%% A is matrix

% clc

A = mprand(4); A = A.'*A;

R = chol(A) %#ok<*NASGU>

L = chol(A,'lower')

R = chol(A,'upper')

try [L,p] = chol(A), catch, fprintf('mp [L,p] = chol(A) BOOOM !!!\n\n'),...

[L,p] = chol(double(A),'lower'), fprintf('matlab ok\n\n'), end %#ok<*ASGLU>

try [L,p] = chol(A,'lower'), catch, fprintf('mp [L,p] = chol(A,''lower'') BOOOM !!!\n\n'),...

[L,p] = chol(double(A),'lower'), fprintf('matlab ok\n\n'), end %#ok<*ASGLU>

try [P,p] = chol(A,'upper'), catch, fprintf('mp [L,p] = chol(A,''upper'') BOOOM !!!\n\n'), ...

[L,p] = chol(double(A),'upper'), fprintf('matlab ok\n\n'), end

%% Why is so important P?

% if your implementation (as I assumed) uses ?potr routine, the value of p will be

% the ?potr's INFO output

% clc

nice = false;

while ~nice

A = rand(4,6);

A = A.'*A;

[~,p] = chol(A);

if p > 2, nice = true; end

end

[L,p] = chol(A)

assert(all(all(L == chol(A(1:p-1,1:p-1)))))

[ma,na] = size(A)

[ma,nl] = size(L)

Please run the following few lines on the last version.

% *** ASSERTIONS *** (This has been incorporated into the documentation in R14SP3):

%

% CHOL expects its input matrix to be symmetric and only looks at the

% upper triangular portion of the matrix.

% CHOL function will return an error if it is only provided with a single

% output argument, and is also given a matrix that is not positive

% definite.

% CHOL function provides an optional second output argument "p" which

% is zero if the matrix is found to be positive definite. If the input

% matrix is not positive definite, then "p" will be a positive integer.

% *** ASSUMPTION ***

% Your implementation uses ?potr routine ...

% >>>>>>>>>> Let's test ...

%% A is scalar

% clc

mp.FollowMatlabNumericFormat(1) % just for convenience

A = mp(9) %#ok<*NOPTS>

R = chol(A) %#ok<*NASGU>

L = chol(A,'lower')

R = chol(A,'upper')

try [L,p] = chol(A,'lower'), catch, fprintf('mp [L,p] = chol(A,''lower'') BOOOM !!!\n\n'),...

[L,p] = chol(double(A),'lower'), fprintf('matlab ok\n\n'), end %#ok<*ASGLU>

try [P,p] = chol(A,'upper'), catch, fprintf('mp [L,p] = chol(A,''upper'') BOOOM !!!\n\n'), ...

[L,p] = chol(double(A),'upper'), fprintf('matlab ok\n\n'), end

%% A is matrix

% clc

A = mprand(4); A = A.'*A;

R = chol(A) %#ok<*NASGU>

L = chol(A,'lower')

R = chol(A,'upper')

try [L,p] = chol(A), catch, fprintf('mp [L,p] = chol(A) BOOOM !!!\n\n'),...

[L,p] = chol(double(A),'lower'), fprintf('matlab ok\n\n'), end %#ok<*ASGLU>

try [L,p] = chol(A,'lower'), catch, fprintf('mp [L,p] = chol(A,''lower'') BOOOM !!!\n\n'),...

[L,p] = chol(double(A),'lower'), fprintf('matlab ok\n\n'), end %#ok<*ASGLU>

try [P,p] = chol(A,'upper'), catch, fprintf('mp [L,p] = chol(A,''upper'') BOOOM !!!\n\n'), ...

[L,p] = chol(double(A),'upper'), fprintf('matlab ok\n\n'), end

%% Why is so important P?

% if your implementation (as I assumed) uses ?potr routine, the value of p will be

% the ?potr's INFO output

% clc

nice = false;

while ~nice

A = rand(4,6);

A = A.'*A;

[~,p] = chol(A);

if p > 2, nice = true; end

end

[L,p] = chol(A)

assert(all(all(L == chol(A(1:p-1,1:p-1)))))

[ma,na] = size(A)

[ma,nl] = size(L)

Kundesupport af UserEcho

***

I have been focusing on covering the MATLAB functionality in breadth, e.g. main functions in prime syntax.

I am always open and happy to do any refinements requested by users.

I believe this user-driven style of development naturally shapes the toolbox functionality, and get audience in different fields.

Besides, that way I am able to actually get some work done:), as some of the functions might take months/years to implement in all possible variations.

I hope for understanding and I am always welcome any requests. Please continue sending your excellent bug reports, comments and suggestions. It is highly appreciated!

***

I plan to release the updated "chol" today or tomorrow the latest....

P.S.

original Mathworks toolboxes rely on p value to check the positiveness definition of matrix.

Here is some example you might find funny (this is pure MATLAB): Matrix A has no positive definite minors, but real positive diagonal tricked MATLAB and it generates unexpected result.

Correct result should be R = [], p =1. Or am I missing something?

The matrix, of course, is not positive defined but the documentation says:

"[R,p] = chol(A)

~~for positive definite A, produces an upper triangular matrix R from the diagonal and upper triangle of matrix A, satisfying the equation R'*R=A and p is zero.~~If A is not positive definite, then p is a positive integer and MATLAB® does not generate an error. When A is full, R is an upper triangular matrix of order q=p-1 such that R'*R=A(1:q,1:q). ..."% make matrix

hermitian

% "... The chol function assumes that A is (complex Hermitian) symmetric. If it is not, chol uses the (complex conjugate) transpose of the upper triangle as the lower triangle ..."

m=triu(A); % [R,p] = chol(A) is the same as [R,p] = chol(A,'upper')

d = diag(m);

m = m + m';

s = size(m,1);

m(1:s+1:s*s) = d;

disp(m)

100 + 0i -1.0151 + 1.4042i -0.12828 + 0.94526i

-1.0151 - 1.4042i 100 + 0i -0.15982 + 1.622i

-0.12828 - 0.94526i -0.15982 - 1.622i -0.12722 + 0.70364i

>> [R,p] = chol(m)

R =

10 + 0i -0.10151 + 0.14042i

0 + 0i 9.9985 + 0i

p =

3

>> R'*R - m(1:p-1,1:p-1)

ans =

0 0

0 1.4211e-14

>> eps(A(p-1,p-1))

ans =

1.4211e-14

****

As implementation I

guessthat the algorithm used is the Bunch-Kaufman factorization of a symmetric matrix by ?sytrf routine ...shorter & better:% make matrix

hermitian% ...

m = triu(A) + triu(A,1)'

But in my example matrix was

notHermitian. Yet MATLAB/zpotrf was not able to recognize this right away and returned p=3, indicating that A(2,2) minor is HPD! Which is not true (matrix A is not even Hermitian!).Moreover returned R gives: Which contradicts the documentation:

When A is full, R is an upper triangular matrix of order q=p-1 such that R'*R=A(1:q,1:q).That was my whole point. I still think correct result should be R=[], p=1, meaning A has no HPD minors :).

What do you think?

(I have traced this behavior down to zpotrf).

example matrix wasYour

notHermitian !(but it does not really matter for matlabchol! :) ... it's the truth )But ( back to Matlab doc):

"

`R = chol(A)`

produces an upper triangular matrix R from the diagonal and upper triangle of matrix A, satisfying the equation R'*R=A.The chol functionc. If it is not, chol uses the (complex conjugate) transpose of the upper triangle as the lower triangle. Matrix A ... "ASSUMESthat A is (complex Hermitian) symmetri(his matrix A)" ...must be positive definite."Matlab chol function NEVER look your ENTIRE matrixbut ONLYtheuppertriangular part(or lower if the second arg is 'lower')of your matrix.

.From this triangular part of YOUR matrix, Matlab build a virtual OWN matrix AND THAT MATRIX IS ALWAYS hermitian(or symmetric if not complex)The R result is always the decomposition ofTHATHERMITIAN matrix.... and, YES, correct result should be R=[], p=1 and (

I checkedyesterday evensymbolic), if you read my post I wrote " The matrix, of course, is not positive defined but the documentation says: ...".More, I suppose, with chances to be fair, the implementation of chol is a package of algorithms that are, when appropriate, improved.

However, the interface (the public contract of chol function) was maintained.

What I say is:

1.

that the result is fair in light of promises made in the matlab documentation and

2.

I never said that sample matrix should bepositive." ...

A is COMPLEX*16 array, dimension (LDA,N)

On entry, the Hermitian matrix A. If UPLO = 'U', the leading

N-by-N upper triangular part of A contains the upper

triangular part of the matrix A,

. If UPLO = 'L', theand the strictly lowertriangular part of A is not referenced

leading N-by-N lower triangular part of A contains the lower

triangular part of the matrix A,

.and the strictly uppertriangular part of A is not referenced

... "

What do you think?

This clears things out. Thank you!

(I used documentation for 2011a, which has a less strict wording - source of my confusion).

(assert( ~isempty ( chol(A, ...))))NOT implies( Aispositive defined )2. but

([~, p] = chol(A, ...) with p > 0)implies( Ais notpositive defined )