How do I create a function script to check the positive definiteness of a a square matrix of any size?

Hey! I'm currently working on a lab where I need to check if a square matrix is positive and definite. My function script needs to

1. accept one sqaure matrix from the calling program,
2. determine if the matrix is positive and definite, and
3. return the result to the calling program - 1 if positive definite, 0 otherwise.

I need to come up with a pattern for the matrices to be check. For a four by four matrix (M) as an example, my code needs to check M(1:1,1:1), M(1:2,1:2), M(1:3,1:3,), M(1:4,1:4), then M(2:2,2:2), M(2,3:2,3), and so on until the last martix M(4:4,4:4). I need help creating a function using for loops! Below is the code I have as of right now, but I don't think it's right. Please help me write some code! Thank you :)

M = [ 1 2 3 4 ; 5 6 7 8 ; 9 10 11 12 ; 13 14 15 16];
x = NAME_lab01(M);
%%
function x = NAME_Lab01(M)
    
    N = length(M)
    
    [r,c] = size (M)
    
    pos_def = 1
    
    for i = 1: N+1-r
        for j = 1:c+r-1
            if det(i,j) > 0
                pos_def = 1
            else
                pos_def = 0
            end
        end
    end
end

[Julia Edwards]

Opps, for the second line, it should read "x = NAME_lab01 (M)"

[Ameer Hamza]

Can you use the built-in functions like eig() or chol()?

https://www.mathworks.com/matlabcentral/answers/101132-how-do-i-determine-if-a-matrix-is-positive-definite-using-matlab

https://www.mathworks.com/help/matlab/math/determine-whether-matrix-is-positive-definite.html

Or do you want to write the algorithm from scratch?

[Julia Edwards]

I am suposed to write the algorithm from scratch. I have changed my code a little, so now it reads as follows:

M = [ 1 2 3 4 ; 5 6 7 8 ; 9 10 11 12 ; 13 14 15 16]

x = NAME_Lab01(M)

%%
function x = NAME_Lab01(M)

N = length(M)

[r,c] = size (M)

pos_def = 1

for i = 1: N+1-r
    for j = i:r
        if det(M(i,j)) > 0
            pos_def = 1
        else
            pos_def = 0
            
        end
    end
end

if pos_def == 1
    disp('1: The matrix is positive definite')
else
    disp('0: The matrix is not positive definite')
end
end

I am now getting an error that says:

Output argument "x" (and maybe others) not assigned during call to "EdwardsJ_Lab1>EdwardsJ_Lab01".

Error in EdwardsJ_Lab1 (line 3)

x = EdwardsJ_Lab01 (M)

What should I do to make my code work?

[per isakson]

Shouldn't it be

function  pos_def  = EdwardsJ_Lab01(M)

?

[Julia Edwards]

I changed my code to that! Now I have code that is running, but it is not working through all the matrices it should be. My code is as follows :

A = [ 1 2 3 4 ; 5 6 7 8 ; 9 10 11 12 ; 13 14 15 16];

B = [ 7 2 ; 1 2];

C = [2 -1 0 ; -1 2 -1 ; 0 -1 2 ]; 

x = EdwardsJ_Lab01 (C);

%%
function [pos_def] = EdwardsJ_Lab01(M)

N = length(M);

[r,c] = size (M);

pos_def = 1;

for i = 1: N+1-r
    for j = i:r
        if det(M(i,j)) > 0
            pos_def = 1
        else
            pos_def = 0
            
        end
    end
end

if pos_def == 1
    fprintf('1: The matrix is positive definite. \n \n ')
else
    fprintf('0: The matrix is not positive definite. \n \n')
end
end

I entered matrix A and B, and got the correct answers, but when I tested matirx C I got the wrong answer. They should all come out to display '1: The matrix is positive definite.' but C outputs '0: The matrix...' I think it is because the for loops are not running through all the determinates. Why is my code not running through all determinates?

[John D'Errico]

Of course, det has nothing to do with telling you if the matrix is positive definite.

For example, the matrix:

A = -eye(2)
A =
    -1     0
     0    -1
     
det(A)
ans =
     1

has a positive determinant, yet is certainly not positive definite, even though it is symmetric.

Perhaps you misunderstand the meaning of positive definite?

[Bruno Luong]

@John I believe IIRC if all principal minors have positive determinants the the matrix is positive.

[Matt J]

Yes, it is Sylvester's criterion,

https://en.wikipedia.org/wiki/Sylvester%27s_criterion

[Bruno Luong]

That's it!

[Julia Edwards]

So is my matrix C actually not positive definite? Meaning my code is correct?

[Matt J]

Your matrix is not positive definite first of all because it is not symmetric. However, your code is by no means correct. det(i,j) makes no sense, for one thing, and throws an error. The det() command should take a single argument (the matrix whos determinant you want to calculate).

[Julia Edwards]

I was given the definition: "A positive definite matrix has several properties, but the one we are interested in today is that the determinants of all the co-diagonal submatrices are positive. A co-diagonal submatrix is a square matrix of any size, contained within the original matrix, that shares the diagonal with the original matrix."

[Bruno Luong]

"Meaning my code is correct?"

Not so fast the Sylvester's criterion only applies for Hermitiam matrix. Your code didn't check that at any moment.

You need to first symmetrize the matrix M = M + M' before checking determinant of principal minors.

You haven't done both.

[Julia Edwards]

I believe that I updated my code to just take one argument, I will attatch it below:

A = [ 1 2 3 4 ; 5 6 7 8 ; 9 10 11 12 ; 13 14 15 16];

B = [ 7 2 ; 1 2];

C = [2 -1 0 ; -1 2 -1 ; 0 -1 2 ]; 

D = [ 1 1 ; 5 1];

x = EdwardsJ_Lab01 (C);

%%
function [pos_def] = EdwardsJ_Lab01(M)

N = length(M);

[r,c] = size (M);

pos_def = 1;

for i = 1: N+1-r
    for j = i:r
        if det(M(i,j)) > 0
            pos_def = 1
        else
            pos_def = 0
            
        end
    end
end

if pos_def == 1
    fprintf('1: The matrix is positive definite. \n \n ')
else
    fprintf('0: The matrix is not positive definite. \n \n')
end
end

[Matt J]

I prefer my one-line solution below, which also checks if the matrix is symmetric.

[Julia Edwards]

Sadly, I do not believe I can use that. I am supposed to write code using for loops and the det() function.

[Matt J]

You do not believe? It sounds like you are not sure. The 3 requirements you posted say nothing about it.

[Bruno Luong]

"A positive definite matrix has several properties, but the one we are interested in today is that the determinants of all the co-diagonal submatrices are positive. A co-diagonal submatrix is a square matrix of any size, contained within the original matrix, that shares the diagonal with the original matrix."

This person who gives you this "defnition" must assume M is symmetric, or makes a mistake. Because if M is not symmetric, it's not true.

The correct definition ths is

M is definite positive if

x'*M*x > 0 for all vector x with |x| > 0.

Now a counter example of what this person wrote:

>> M=[10 20; 1 10] % not symmetric

M =

    10    20
     1    10

>> det(M(1,1))

ans =

    10

>> det(M(1:2,1:2))

ans =

    80

So is M positive according to this? The answer is no

>> v=[1; -1]

v =

     1
    -1

>> v'*M*v

ans =

    -1 % NEGATIVE

[Julia Edwards]

Why isn't M positive according to your example? Both determinates are positive? The matrix is also symmetrical.

[Bruno Luong]

'The matrix is also symmetrical' : No you should read again my example.

M is not definite postive because v'*M*v = -1 < 0, that violates the definition here

[Paul]

Matt J,

“Your matrix is not positive definite first of all because it is not symmetric.”

I don’t think M symmetric is a necessary criterion for M to be positive definite.

[Bruno Luong]

+1

In this page from Wolfram where I post the link to the definition they wrote:

"Confusingly, the discussion of positive definite matrices is often restricted to only Hermitian matrices, or symmetric matrices in the case of real matrices"

This thread proves they are right again.

[Paul]

Ironically, that same MathWorld page also says

“The definition of positive definiteness is equivalent to the requirement that the determinants associated with all upper-left submatrices are positive.“

But apparently this is true only for symmetric, positive definite matrices.

[Matt J]

I have never seen a text that does not predicate the definition of positive definiteness on (Hermitian) symmetry,

[Bruno Luong]

I'll show you one more:

Here is a definition from the famous MATHWORKS support team, "A non-symmetric matrix (B) is positive definite if all eigenvalues of (B+B')/2 are positive."

PS: John Petersen gives there also a count example of non-symmetric non-DP matrix yet all principal minors has positive determinant.

[Matt J]

With all due respect to TMW, I don't know if I would equate their definitions with textbook definitions.

[Paul]

I don’t have any texts nearby but will check some when I can out of curiosity. The Wikipedia page does talk about the non-Hermitian case in the “Extension ...” section. Interestingly, that Wikipedia page has an external link to the MathWorld page that uses the more general definition.

My recollection is that because the quadratic form of M only depends on the symmetric part of M, it suffices to assume M is symmetric, because if it isn’t one can always use just the symmetric part of M for analysis of its quadratic form anyway.

[Bruno Luong]

Golub Van Loan definition, you know them?

TMW (Matlab), Wolfram (Mathematica), Golub/Van-Loan (the bible matrix computations). It seems enough authorities to me concerning linear algebra.

[Paul]

MathWorld and GVL are sufficient for me to stick with the less restrictive definition that I use.

[Matt J]

Ah well, if Golub and Van-Loan says so, then fine. I've never read the bible. I'm not that religious... But do non-symmetric positive definite matrices have any useful properties? it seems like you need the symmetry assumption to apply Sylvester's criterion, the eigenvalue criterion, and pretty much every other property I've seen used in practice.

[Paul]

I don’t know if NSPD matrices have additional useful properties compared to SPD matrixes. But they must have some importance in solving linear systems as in the GVL link that Bruno posted.

Maybe I’m just parsing words here, but you don’t need the symmetry assumption to apply those criteria. Rather the PD of M, symmetric or not, can always be determined by applying them to M+M’ . More generally, I’m pretty sure that everything you need to know about the quadratic form of M can be determined from the quadratic form of (M+M’)/2.

[Bruno Luong]

Personally I never deal with unsymmetric DP matrix, but when I was taugh bilinear form, they teach us a bilinear for can be NOT necessary symmetric, which represents by an unsymmetric matrix.

I guess the theory can also be extrended to some binear form that is applied on vector of elements that belong non-cummutative ring (such as quaternion), in which you can't swap the order.

Quantum mechanics manipulates typical such bilinear form, but I'm unable to tell you a concrete example of operator that is unsymmetric and DP, this is not my field of expertise.

But then Golub/Van-Loan showed in this paper they can exploit such definite positive property on non-symeetric matrix and doing some improvement on a basic task of solving the associate linear equations.

Answer 1

[Matt J]

If you're allowed to use stock Matlab commands like det(), I don't know why you couldn't use eig().

pos_def = isequal(M,M') && all( eig(M)>0 )

Answer 2

[Bruno Luong]

Code based on Sylvester's critterion

function tf = isposdef(M)
tf = ispd(M+M');
end

function tf = ispd(M)
tf = det(M)>0 && (length(M)<=1 || ispd(M(1:end-1,1:end-1)));
end

EDIT: fix bug length(M)<=1

[Matt J]

The logic in ispd isn't right. The code incorrectly tells me that

M=[-4 1;1 -1];

is positive definite.

[Matt J]

I think it should be

tf = det(M)>0 && (length(M)>=1 && ispd(M(1:end-1,1:end-1)));

[Bruno Luong]

Yes I edited by hand and made a mistake. Fix now.

[Paul]

Bruno,

I'm sure you know this, but I'm going to post this comment for potential interest of others.

The input argument M to ispd is Hermitian and therefore theoretically its determinant is real. However, due to roundoff errors the computed value of det(M) could have a small imaginary part, which got me thinking about what det(M) > 0 means if det(M) returns a complex number. Sure enough, the gt operator only considers the real part of a complex operand.

[Bruno Luong]

Paul, yes very good comment. I'm awared of MATLAB z > 0 is equivalent to real(z) > 0 for z complex. So the code is robust to numerical error without the need of cheking real(det(M)) > 0.

Though the best method to check PD IMO is using incomplete cholesky factorization for full and EIGS for sparse matrix, and not using DET at all.

[Paul]

Is the incomplete Cholesky factorization implemented by the chol function? I'm asking because of the thread in the first link that Ameer posted way earlier in this thread that shows that the chol function may not be reliable in determining if a matrix is SPD (repeating the link here for easy access): https://www.mathworks.com/matlabcentral/answers/101132-how-do-i-determine-if-a-matrix-is-positive-definite-using-matlab?s_tid=srchtitle

Repeating an example from that thread:

A =0.0327*ones(2,2);
[R,flag]=chol(A);flag

flag =

     0
     

According to doc chol (R2019A): "If flag = 0 then the input matrix is symmetric positive definite and the factorization was successful." In this case, A is definitely symmetric, but I don't think it's PD.

[Bruno Luong]

IMO there is no algorithm that can 100% reliable in case the matrix is non-negative and not definite positive (one eigen value is 0), since any method will be sensitive to numerical errors that can randomly make or break the result.

But yeah intuitively I would also think CHOL could be less reliable than EIGS.

Answer 3

[Matt J]

Sadly, I do not believe I can use that. I am supposed to write code using for loops and the det() function.

The code below fulfills the "for-loops +det() function" requirement. This is to make the point that the assignment is non-sensical unles it specifically says you must use Sylvester's criterion.

function tf=isposdef(M)

     e=eig(M+M');
     
     tf=true;
     for i=1:numel(e)
         tf=tf & det(e(i))>0;
     end

end