0

Completat

chol (~call) again

fa 9 anys • actualitzat fa 9 anys • 15

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)

Vote

Replies 15
vells primer
- nous primer
- vells primer

Started

fa 9 anys

Yes, up to this moment we supported only the following syntax:

R = chol(A)
L = chol(A,'lower')
R = chol(A,'upper')

I understand that the two-outputs version is important and it will be ready shortly :).

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

Respondre
|

fa 9 anys

As for me, no hurry! my mp.m is already patched :)

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

Respondre
|

fa 9 anys

Just completed adding the feature - not building the release.

Here is some example you might find funny (this is pure MATLAB):

>> A =
          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(A)
R =
           10 +          0i     -0.10151 +    0.14042i
            0 +          0i       9.9985 +          0i
p =
     3

>> norm(R'*R-A(1:p-1,1:p-1),1)
ans =
       2.8084

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?

Respondre
|

fa 9 anys

interesting! I like! I'll think about it ...

Respondre
|

fa 9 anys

No, Matlab was not tricked; I'm disappointed :); The result is according to the documentation!

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 guess that the algorithm used is the Bunch-Kaufman factorization of a symmetric matrix by ?sytrf routine ...

Respondre
|

fa 9 anys

shorter & better:
% make matrix hermitian
% ...
m = triu(A) + triu(A,1)'

Respondre
|

fa 9 anys

Yes, your matrix is Hermitian, of course.

But in my example matrix was not Hermitian. 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:

>> norm(R'*R-A(1:p-1,1:p-1),1)
ans =
       2.8084

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

Respondre
|

fa 9 anys

Motto: "No matter what you want, what you get counts"

Your example matrix was not Hermitian ! (but it does not really matter for matlab chol! :) ... 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 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. Matrix A ... " (his matrix A) " ...must be positive definite."

Matlab chol function NEVER look your ENTIRE matrix but ONLY the upper triangular 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 of THAT HERMITIAN matrix.

... and, YES, correct result should be R=[], p=1 and (I checked yesterday even symbolic), 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 be positive.

Respondre
|

fa 9 anys

I guess I'm missing something ... what does the norm matrix with problems raised by you, right or wrong behavior of the Matlab function chol?

Respondre
|

fa 9 anys

From zpotrf.f

" ...
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, and the strictly lower
triangular part of A is not referenced. If UPLO = 'L', the
leading N-by-N lower triangular part of A contains the lower
triangular part of the matrix A, and the strictly upper
triangular part of A is not referenced.
... "

What do you think?

Respondre
|

fa 9 anys

I think you are right. For some reason I assumed that second output from chol should also reflect if matrix Hermitian or not. But it does check only positive definiteness.

This clears things out. Thank you!
(I used documentation for 2011a, which has a less strict wording - source of my confusion).

Respondre
|

fa 9 anys

if I had a thousand dollars for each data when I'm wrong, that would be great :) But how I'm not Microsoft employee's, I do not receive the money.

Respondre
|

fa 9 anys

1.
(assert( ~isempty ( chol(A, ...)))) NOT implies ( A is positive defined )
2. but
([~, p] = chol(A, ...) with p > 0) implies ( A is not positive defined )

Respondre
|

fa 9 anys

New version has been released, please download it from our website.

Respondre
|

Completat

fa 9 anys

Respondre
|

Customer support service by UserEcho