1. General
  2. Ideed
0
Lõpetatud

chol (~call) again

john doe 9 aastat tagasi • uuendatud 9 aastat tagasi • 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)







Alustatud
Pavel Holoborodko 9 aastat tagasi
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....
john doe 9 aastat tagasi
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.
Pavel Holoborodko 9 aastat tagasi
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?
john doe 9 aastat tagasi
interesting! I like! I'll think about it ...
john doe 9 aastat tagasi
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 ...
john doe 9 aastat tagasi
shorter & better:
% make matrix hermitian
% ...
m = triu(A) + triu(A,1)'
Pavel Holoborodko 9 aastat tagasi
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).

john doe 9 aastat tagasi
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.
john doe 9 aastat tagasi
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?
john doe 9 aastat tagasi
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?
Pavel Holoborodko 9 aastat tagasi
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).

john doe 9 aastat tagasi
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.
john doe 9 aastat tagasi
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 )
Pavel Holoborodko 9 aastat tagasi
New version has been released, please download it from our website.
Lõpetatud
Pavel Holoborodko 9 aastat tagasi

Customer support service by UserEcho