Is the Matrix Positive-Definite?

19

Introduction

Today we're gonna take care of the bane of first-year linear algebra students: matrix definiteness! Apparently this doesn't yet have a challenge so here we go:

Input

  • A \$n\times n\$ symmetric Matrix \$A\$ in any convenient format (you may also of course only take the upper or the lower part of the matrix)
  • Optionally: the size of the matrix \$n\$

What to do?

The challenge is simple: Given a real-valued matrix \$n\times n\$ Matrix decide whether it is positive definite by outputting a truthy value if so and a falsey value if not.

You may assume your built-ins to actually work precisely and thus don't have to account for numerical issues which could lead to the wrong behaviour if the strategy / code "provably" should yield the correct result.

Who wins?

This is , so the shortest code in bytes (per-language) wins!


What is a positive-definite Matrix anyways?

There are apparently 6 equivalent formulations of when a symmetric matrix is positive-definite. I shall reproduce the three easier ones and reference you to Wikipedia for the more complex ones.

  • If \$\forall v\in\mathbb R^n\setminus \{0\}: v^T Av>0\$ then \$A\$ is positive-definite.
    This can be re-formulated as:
    If for every non-zero vector \$v\$ the (standard) dot product of \$v\$ and \$Av\$ is positive then \$A\$ is positive-definite.
  • Let \$\lambda_i\quad i\in\{1,\ldots,n\}\$ be the eigenvalues of \$A\$, if now \$\forall i\in\{1,\ldots,n\}:\lambda_i>0\$ (that is all eigenvalues are positive) then \$A\$ is positive-definite.
    If you don't know what eigenvalues are I suggest you use your favourite search engine to find out, because the explanation (and the needed computation strategies) is too long to be contained in this post.
  • If the Cholesky-Decomposition of \$A\$ exists, i.e. there exists a lower-triangular matrix \$L\$ such that \$LL^T=A\$ then \$A\$ is positive-definite. Note that this is equivalent to early-returning "false" if at any point the computation of the root during the algorithm fails due to a negative argument.

Examples

For truthy output

\begin{pmatrix}1&0&0\\0&1&0\\0&0&1\end{pmatrix}

\begin{pmatrix}1&0&0&0\\0&2&0&0\\0&0&3&0\\0&0&0&4\end{pmatrix}

\begin{pmatrix}5&2&-1\\2&1&-1\\-1&-1&3\end{pmatrix}

\begin{pmatrix}1&-2&2\\-2&5&0\\2&0&30\end{pmatrix}

\begin{pmatrix}7.15&2.45\\2.45&9.37\end{pmatrix}

For falsey output

(at least one eigenvalue is 0 / positive semi-definite) \begin{pmatrix}3&-2&2\\-2&4&0\\2&0&2\end{pmatrix}

(eigenvalues have different signs / indefinite) \begin{pmatrix}1&0&0\\0&-1&0\\0&0&1\end{pmatrix}

(all eigenvalues smaller than 0 / negative definite) \begin{pmatrix}-1&0&0\\0&-1&0\\0&0&-1\end{pmatrix}

(all eigenvalues smaller than 0 / negative definite) \begin{pmatrix}-2&3&0\\3&-5&0\\0&0&-1\end{pmatrix}

(all eigenvalues smaller than 0 / negative definite) \begin{pmatrix}-7.15&-2.45\\-2.45&-9.37\end{pmatrix}

(three positive, one negative eigenvalue / indefinite) \begin{pmatrix}7.15&2.45&1.23&3.5\\2.45&9.37&2.71&3.14\\1.23&2.71&0&6.2\\3.5&3.14&6.2&0.56\end{pmatrix}

SEJPM

Posted 2018-09-23T11:47:44.397

Reputation: 3 203

This challenge was sandboxed. – SEJPM – 2018-09-23T11:48:06.563

You need to provide a better definition of what we're supposed to be looking for rather than assuming we can all read mathematical notation (or all know what an "eigenvalue" is). A worked example would be useful too. – Shaggy – 2018-09-23T15:53:04.303

9@Shaggy I think the challenge is better without all the background to clog it up. There are many existing explanations of what an eigenvalue is elsewhere, and this post is already really large. – Post Rock Garf Hunter – 2018-09-23T16:37:11.803

@Shaggy I won't add a self-contained description / introduction of eigenvalues because this is neither the right place for that nor am I probably qualified enough to do that in a sufficient manner. I will also not add a "worked example" because there is no single "superior" strategy to this problem, because right now answers are generally using three different approaches to this challenge (either finding the eigenvalues and then checking they're positive or checking whether the cholesky decomposition exists or Sylvester's criterion). – SEJPM – 2018-09-23T17:40:26.320

1The challenge would have been nicer hadn't you restrict the input to symmetric matrices. – polfosol ఠ_ఠ – 2018-09-24T15:55:00.007

@SEJPM yes I noticed it. That was the point of my comment – polfosol ఠ_ఠ – 2018-09-24T16:07:03.693

@polfosolఠ_ఠ oh, I misread your comment to be a complaint that symmetric matrices are not guaranteed. Also I found it to be a boring tag-on to also force people to check for symmetric-ness. – SEJPM – 2018-09-24T16:09:50.130

1I meant just checking for the sign of eigenvalues is also boring. Different tastes I know ;) – polfosol ఠ_ఠ – 2018-09-24T16:11:10.073

Answers

11

C, 108 bytes

-1 byte thanks to Logern
-3 bytes thanks to ceilingcat

f(M,n,i)double**M;{for(i=n*n;i--;)M[i/n][i%n]-=M[n][i%n]*M[i/n][n]/M[n][n];return M[n][n]>0&(!n||f(M,n-1));}

Try it online!

Performs Gaussian elimination and checks whether all diagonal elements are positive (Sylvester's criterion). Argument n is the size of the matrix minus one.

nwellnhof

Posted 2018-09-23T11:47:44.397

Reputation: 10 037

Perhaps save a character with float instead of double? – Jens – 2018-09-23T15:10:30.950

106 bytes - Try it online! – Logern – 2018-09-23T15:35:05.423

You can shave another character if you drop i=0 in the for loop, make the recursive call f(M,n-1,0) and the initial call with 0 as the third arg. – Jens – 2018-09-23T19:36:17.190

@Jens 1. Using floats instead of doubles can quickly lead to noticable rounding errors, so I don't think the one byte saved is worth it. 2. Initializing a variable via an additional argument looks like cheating to me. – nwellnhof – 2018-09-24T17:42:33.150

@Logern I refuse to use the "omit the return statement" trick in my C answers. But thanks for the other byte saved. – nwellnhof – 2018-09-24T17:45:56.837

107 bytes - Try it online! – Logern – 2018-09-24T18:18:13.147

9

MATLAB/Octave, 19 17 12 bytes

@(A)eig(A)>0

Try it online!

The function eig provides the eigenvalues in ascending order, so if the first eigenvalue is greater than zero, the other ones are too.

Daniel Turizo

Posted 2018-09-23T11:47:44.397

Reputation: 101

You can drop the f= at the beginning - anonymous functions are generally accepted as answers. – Delfad0r – 2018-09-23T15:05:21.880

Thanks for the tip! – Daniel Turizo – 2018-09-23T15:11:53.473

Even if its a vector? Interesting – Daniel Turizo – 2018-09-23T15:17:55.363

https://codegolf.meta.stackexchange.com/a/2194/82619 – Delfad0r – 2018-09-23T15:21:46.193

1+1. I've added a link for trying it online. Hope you don't mind. Note that it is also proving that the output values despite being arrays do count as the correct "truthy" or "falsey" values as per the link @Delfad0r posted. – Tom Carpenter – 2018-09-23T17:03:22.330

2Having said that, it fails for the first "falsey" test case on TIO. I'm guessing due to a precision issue - one of the Eigen values comes out as 8.9219e-17 rather than 0. – Tom Carpenter – 2018-09-23T17:05:43.017

Maybe eig(a)>eps could help. But I don't think its a good idea. – tsh – 2018-09-23T17:10:12.370

@TomCarpenter thank you! It is in fact a precision issue: if you input the matrix as symbolic, it will yield the correct result. Altough it must be noticed that for some sym matrices the result is not a bool vector but a sym vector. – Daniel Turizo – 2018-09-23T17:32:19.470

Not that your answer depends on it, but the Matlab help page for "eig" says that it doesn't guarantee the order which eigenvalues appear (certainly for non-symmetric matrices, you'll get them in non-sorted order; the algorithm currently used for symmetric matrices when detected may impose sorting). – Batman – 2018-09-25T02:41:08.967

7

Jelly, 11 10 bytes

ṖṖ€$ƬÆḊṂ>0

Uses Sylvester's criterion.

Try it online!

How it works

ṖṖ€$ƬÆḊṂ>0  Main link. Argument: M (matrix)

   $Ƭ       Do the following until a fixed point is encountered.
Ṗ             Pop; remove the last row of the matrix.
 Ṗ€           Pop each; remove the last entry of each row.
     ÆḊ     Take the determinants of the resulting minors.
       Ṃ    Take the minimum.
        >0  Test if the least determinant is positive, i.e., if all determinants are.

Dennis

Posted 2018-09-23T11:47:44.397

Reputation: 196 637

6

R, 29 bytes

function(m)all(eigen(m)$va>0)

Try it online!


Alternative using cholesky :

R, 34 33 bytes

function(m)is.array(try(chol(m)))

Try it online!

-1 byte thanks to @Giuseppe

digEmAll

Posted 2018-09-23T11:47:44.397

Reputation: 4 599

6

Haskell, 56 bytes

f((x:y):z)=x>0&&f[zipWith(-)v$map(u/x*)y|u:v<-z]
f[]=1>0

Try it online!

Basically a port of nwellnhof's answer. Performs gaussian elimination and checks whether the elements on the main diagonal are positive.

Fails the first falsey output because of rounding errors, but it would theoretically work with infinite precision. Thanks to Curtis Bechtel's suggestion, now the outputs are all correct.

Delfad0r

Posted 2018-09-23T11:47:44.397

Reputation: 1 688

2you can add inputs :: [[[Rational]]] to get correct answers – Curtis Bechtel – 2018-09-23T20:38:46.583

4

Wolfram Language (Mathematica), 20 bytes

0<Min@Eigenvalues@#&

Try it online!

Mr. Xcoder

Posted 2018-09-23T11:47:44.397

Reputation: 39 774

Should 4th test case be False? – tsh – 2018-09-23T14:48:47.030

@tsh Fixed, I am dumb! – Mr. Xcoder – 2018-09-23T14:57:40.890

8

Funny how Mathematica has a builtin for this, but its name is longer than your solution.

– Federico Poloni – 2018-09-23T15:04:36.300

@FedericoPoloni: Wouldn't a solution using NullSpace or MatrixRank be even shorter? If Null space is zero then the matrix is positive definite. – Phil H – 2018-09-24T15:41:47.697

@PhilH No, I'm afraid that doesn't work by itself. For instance, the second falsey example (diagonal matrix with (1,-1,1)) has rank 3, but is not positive definite. – Federico Poloni – 2018-09-24T16:39:21.007

3

MATL, 4 bytes

Yv0>

Try it online!

I've noticed that for the [3 -2 2; -2 4 0; 2 0 2] test case, MATL calculates the 0 eigenvalue with a very small inaccuracy, computing something as low as about \$10^{-18}\$. This therefore fails the first falsy test case due to precision issues. Thanks to Luis Mendo for pointing out that a non-empty array is truthy iff all its entries differ from 0, saving 1 byte.

Mr. Xcoder

Posted 2018-09-23T11:47:44.397

Reputation: 39 774

1@LuisMendo Thanks, I learned something new about MATL today! – Mr. Xcoder – 2018-09-26T21:53:46.423

My pleasure :-) Here's a better explanation by Suever. I forgot to say that for complex-valued arrays only the real part is compared against zero. So [1 2 3j] is falsey

– Luis Mendo – 2018-09-26T22:04:48.963

2

Mathematica, 29 28 bytes

AllTrue[Eigenvalues@#,#>0&]&

Definition 2.

Shieru Asakoto

Posted 2018-09-23T11:47:44.397

Reputation: 4 445

2

MATL, 6 bytes

It is possible to do it using even less bytes, @Mr. Xcoder managed to find a 5 byte MATL answer!

YvX<0>

Explanation

Yv     compute eigenvalues
  X<   take the minimum
    0> check whether it is greather than zero

Try it online!

flawr

Posted 2018-09-23T11:47:44.397

Reputation: 40 560

Fails on the first falsy test case. See my deleted answer. – Mr. Xcoder – 2018-09-25T06:59:48.953

1@Mr.Xcoder Oh, you even submitted an answer before me. I think you should undelete your answer as it just depends on rounding issues. (I think you can expect answers to use limited precision arithmetic - I think only the CAS languages here use exact computations.) – flawr – 2018-09-25T18:02:33.147

Following your advice, I've undeleted it.

– Mr. Xcoder – 2018-09-25T21:18:54.977

2

Julia 1.0, 28 bytes

using LinearAlgebra
isposdef

Try it online!


Julia 0.6, 8 bytes

isposdef

Try it online!

Lyndon White

Posted 2018-09-23T11:47:44.397

Reputation: 1 021

1

Maple, 33 bytes

(i.e. my 2 cents)

with(LinearAlgebra):
IsDefinite(A)

polfosol ఠ_ఠ

Posted 2018-09-23T11:47:44.397

Reputation: 519

Hello and welcome to PPCG; I am unfamiliar with Maple, though is the newline necessary? – Jonathan Frech – 2018-09-24T16:06:12.440

@JonathanFrech Hello and thanks. No it's not. I didn't count it btw. – polfosol ఠ_ఠ – 2018-09-24T16:09:06.233

To me your current byte count seems to reflect the newline character. – Jonathan Frech – 2018-09-24T16:16:17.700

@JonathanFrech Duh, my bad – polfosol ఠ_ఠ – 2018-09-24T16:18:13.160

1Well ... now your code and your byte count disagree. – Jonathan Frech – 2018-09-24T17:07:36.890

0

JavaScript (ES6),  99 95  88 bytes

Same method as nwellnhof's answer. Returns \$0\$ or \$1\$.

f=(m,n=0,R=m[n])=>R?f(m,n+1)&R[m.map((r,y)=>y<n&&R.map((v,x)=>r[x]-=v*r[n]/R[n])),n]>0:1

Try it online!

Arnauld

Posted 2018-09-23T11:47:44.397

Reputation: 111 334