I know the definition of symmetric positive definite (SPD) matrix, but want to understand more.

Why are they so important, intuitively?

Here is what I know. What else?

  • For a given data, Co-variance matrix is SPD. Co-variance matrix is a important metric, see this excellent post for intuitive explanation.

  • The quadratic form $\frac 1 2 x^\top Ax-b^\top x +c$ is convex, if $A$ is SPD. Convexity is a nice property for a function that can make sure the local solution is global solution. For Convex problems, there are many good algorithms to solve, but not for non-covex problems.

  • When $A$ is SPD, the optimization solution for the quadratic form $$\text{minimize}~~~ \frac 1 2 x^\top Ax-b^\top x +c$$ and the solution for linear system $$Ax=b$$ are the same. So we can run conversions between two classical problems. This is important because it enables us to use tricks discovered in one domain in the another. For example, we can use the conjugate gradient method to solve a linear system.

  • There are many good algorithms (fast, numerical stable) that work better for an SPD matrix, such as Cholesky decomposition.

EDIT: I am not trying ask the identities for SPD matrix, but the intuition behind the property to show the importance. For example, as mentioned by @Matthew Drury, if a matrix is SPD, Eigenvalues are all positive real numbers, but why all positive matters. @Matthew Drury had a great answer to flow and that is what I was looking for.

  • 5
    Eigenvalues are all positive real numbers. This fact underlies many of the others. – Matthew Drury Jul 15 '16 at 19:22
  • 3
    To go a little further than @Matthew: If you choose a suitable basis, all such matrices are the same and are equal to the identity matrix. In other words, there is exactly one positive-definite quadratic form in each dimension (for Real vector spaces) and it's the same as Euclidean distance. – whuber Jul 15 '16 at 19:55
  • 2
    You'll find some intuition in the many elementary ways of showing the eigenvalues of a real symmetric matrix are all real: mathoverflow.net/questions/118626/… In particular, the quadratic form $x^TAx$ occurs naturally in the Rayleigh quotient, and symmetric matrices provide natural way of exhibiting a large family of matrices whose eigenvalues are real. See the Courant minimax theorem for example: en.wikipedia.org/wiki/Courant_minimax_principle – Alex R. Jul 15 '16 at 23:28
  • 4
    This seems overly broad;if it didn't already have three answers I'd likely have closed it on that basis. Please offer more guidance about what you specifically want to know (asking for intuition is much too personal/individual for people to guess at in a case like this) – Glen_b Jul 16 '16 at 4:05
  • 1
    I am having a hard time of coming up a situation in statistics that would give rise to a matrix that is not p.s.d. (unless you screwed up in computing a correlation matrix, e.g. by filling it up with pairwise correlation computed on data with missing values). Any square symmetric matrix I can think of is either a covariance, an information or a projection matrix. (Elsewhere in applied mathematics, the non-p.s.d. matrices may be a cultural norm, e.g. the finite element matrices in PDE, say.) – StasK Aug 11 '16 at 17:08
up vote 13 down vote accepted

A (real) symmetric matrix has a complete set of orthogonal eigenvectors for which the corresponding eigenvalues are are all real numbers. For non-symmetric matrices this can fail. For example, a rotation in two dimensional space has no eigenvector or eigenvalues in the real numbers, you must pass to a vector space over the complex numbers to find them.

If the matrix is additionally positive definite, then these eigenvalues are all positive real numbers. This fact is much easier than the first, for if $v$ is an eigenvector with unit length, and $\lambda$ the corresponding eigenvalue, then

$$ \lambda = \lambda v^t v = v^t A v > 0 $$

where the last equality uses the definition of positive definiteness.

The importance here for intuition is that the eigenvectors and eigenvalues of a linear transformation describe the coordinate system in which the transformation is most easily understood. A linear transformation can be very difficult to understand in a "natural" basis like the standard coordinate system, but each comes with a "preferred" basis of eigenvectors in which the transformation acts as a scaling in all directions. This makes the geometry of the transformation much easier to understand.

For example, the second derivative test for the local extrema of a function $R^2 \rightarrow R$ is often given as a series of mysterious conditions involving an entry in the second derivative matrix and some determinants. In fact, these conditions simply encode the following geometric observation:

  • If the matrix of second derivatives is positive definite, you're at a local minimum.
  • If the matrix of second derivatives is negative definite, you're at a local maximum.
  • Otherwise, you are at neither, a saddle point.

You can understand this with the geometric reasoning above in an eigenbasis. The first derivative at a critical point vanishes, so the rates of change of the function here are controlled by the second derivative. Now we can reason geometrically

  • In the first case there are two eigen-directions, and if you move along either the function increases.
  • In the second, two eigen-directions, and if you move in either the function decreases.
  • In the last, there are two eigen-directions, but in one of them the function increases, and in the other it decreases.

Since the eigenvectors span the whole space, any other direction is a linear combination of eigen-directions, so the rates of change in those directions are linear combinations of the rates of change in the eigen directions. So in fact, this holds in all directions (this is more or less what it means for a function defined on a higher dimensional space to be differentiable). Now if you draw a little picture in your head, this makes a lot of sense out of something that is quite mysterious in beginner calculus texts.

This applies directly to one of your bullet points

The quadratic form $\frac 1 2 x^\top Ax-b^\top x +c$ is convex, if $A$ is SPD. Convex is a nice property that can make sure the local solution is global solution

The matrix of second derivatives is $A$ everywhere, which is symmetric positive definite. Geometrically, this means that if we move away in any eigen-direction (and hence any direction, because any other is a linear combination of eigen-directions) the function itself will bend away above it's tangent plane. This means the whole surface is convex.

  • 5
    A graphical way of looking at it: if $\mathbf A$ is SPD, the contours of the associated quadratic form are ellipsoidal. – J. M. is not a statistician Jul 16 '16 at 14:06
  • 7
    That characterization by @J.M. is very perceptive. In case anyone is wondering what might be special about ellipsoidal contours, note that they are just perfect spheres in disguise: the units of measurement may differ along their principal axes and the ellipsoids might be rotated with respect to the coordinates in which the data are described, but for a great many purposes--especially conceptual ones--those differences are inconsequential. – whuber Jul 18 '16 at 16:06
  • That's related to my way of understanding Newton's method geometrically. Best approximate the current level set with an ellipsoid, and then take a coordinate system where the ellipsoid is a circle, move orthogonal to the circle in that coordinate system. – Matthew Drury Jul 18 '16 at 16:13
  • @whuber And if we want to make that general assessment more specific, one could point out that this is what happens when we rotate data into a PCA coordinate system -- among many other types of transformations. (Do I have that right?) – Sycorax Jul 18 '16 at 17:36
  • 1
    If there are (active) constraints, you need to project into the Jacobian of the active constraints before doing the eigenvalue and eigendirection spiel. If the Hessian is psd, the (any) projection will be psd, but the converse is not necessarily true, and often isn't. See my answer. – Mark L. Stone Oct 2 '16 at 14:00

You'll find some intuition in the many elementary ways of showing the eigenvalues of a real symmetric matrix are all real: https://mathoverflow.net/questions/118626/real-symmetric-matrix-has-real-eigenvalues-elementary-proof/118640#118640

In particular, the quadratic form $x^TAx$ occurs naturally in the Rayleigh quotient, and symmetric matrices provide what's arguably the most natural way of exhibiting a large family of matrices whose eigenvalues are real. See the Courant minimax theorem for example: https://en.wikipedia.org/wiki/Courant_minimax_principle

Also symmetric, strictly positive definite matrices are the only set of matrices which can define a non-trivial inner product, along with an induced norm: $d(x,y)=\langle x,Ay\rangle=x^TAy$. This is because by definition for real vectors $x,y$ $d(x,y)=d(y,x)$ for all $x,y$ and $\|x\|^2=x^TAx>0$ for $x\neq 0$. In this way, symmetric positive definite matrices can be viewed as ideal candidates for coordinate transforms.

This latter property is absolutely key in the area of support vector machines , specifically kernel methods and the kernel trick, where the kernel must be symmetric positive to induce the right inner product. Indeed Mercer's theorem generalizes the intuitive properties of symmetric matrices to functional spaces.

Geometrically, a positive definite matrix defines a metric, for instance a Riemannian metric, so we can immediately use geometric concepts.

If $x$ and $y$ are vectors and $A$ a positive definite matrix, then $$ d(x,y) = (x-y)^T A (x-y) $$ is a metric (also called distance function).

In addition, positive definite matrices are related to iner product: In $\mathbb{R}^n$, we can define an inner product by $$ \langle x,y \rangle = x^T A y $$ where $A$ as above is positive definite. More, all inner products on $\mathbb{R}^n$ arises in this way.

With respect to optimization (because you tagged your question with the optimization tag), SPD matrices are extremely important for one simple reason - an SPD Hessian guarantees that the search direction is a descent direction. Consider the derivation of Newton's method for unconstrained optimization. First, we form the Taylor expansion of $f(x + \Delta x)$:

$$f(x + \Delta x)\approx f(x) + \Delta x^T \nabla f(x)+ \frac{1}{2} \Delta x^T \nabla^2 f(x) \Delta x$$

Next, we take the derivative with respect to $\Delta x$:

$$f'(x + \Delta x)\approx \nabla f(x) + \nabla^2 f(x) \Delta x$$

Finally, set the derivative equal to 0 and solve for $\Delta x$:

$$\Delta x = -\nabla^2 f(x)^{-1} \nabla f(x)$$

Assuming $\nabla^2 f(x)$ is SPD, it is easy to see that $\Delta x$ is a descent direction because:

$$\nabla f(x)^T \Delta x = -\nabla f(x)^T \nabla^2 f(x)^{-1} \nabla f(x) < 0$$

When using Newton's method, non-SPD Hessian matrices are typically "nudged" to be SPD. There's a neat algorithm called modified Cholesky that will detect a non-SPD Hessian, "nudge" it appropriately in the right direction and factorize the result, all for (essentially) the same cost as a Cholesky factorization. Quasi-Newton methods avoid this problem by forcing the approximate Hessian to be SPD.

As an aside, symmetric indefinite systems are receiving a lot of attention these days. They come up in the context of interior point methods for constrained optimization.

  • Thank you very much for great answer. I understand decent direction is important in line search method. In trust region methods, decent direction is also important? – hxd1011 Jul 18 '16 at 17:06
  • 1
    It is still important for trust region methods. Trust region methods basically work by bounding the step size FIRST and then solving for the step direction. If the step does not achieve the desired decrease in objective function value, you reduce the bound on the step size and start over. Imagine that your algorithm for generating the step direction does not guarantee that the step direction is a descent direction. Even as the radius of the trust region goes to 0, you may never generate an acceptable step (even if one exists) because none of your step directions are descent directions. – Bill Woessner Jul 19 '16 at 17:36
  • Line search methods basically exhibit the same behavior. If your search direction is not a descent direction, the line search algorithm may never find an acceptable step length - because there isn't one. :-) – Bill Woessner Jul 19 '16 at 17:38
  • Great answer, thank you for helping me to connect the pieces. – hxd1011 Jul 19 '16 at 17:40

An often under-appreciated trick is that symmetric matrices easily permit one to compute the matrix exponential. While it is true in general that there are "19 Dubious Ways to Compute the Matrix Exponential" (Cleve Moler & Charles Van Loan), symmetric matrices, and therefore SPD matrices, have the property that the spectral decomposition can be used to compute the matrix exponential precisely and in a numerically-stable fashion.* Consider the factorization $A=PDP^{-1}$, where $P$ is a matrix of eigenvectors and $D$ is a diagonal matrix of eigenvalues. Because $A$ is symmetric, that means $P$ is orthonormal: $P^T=P^{-1}$. The matrix exponential may be computed as $\exp(A)=P\exp(D)P^T,$ where $D$ may be cheaply computed by exponentiating all elements along the diagonal. Proofs of this method can be convstructed by the definition of the exponential function as a power series.

Indeed, many other computations can be computed in the same way. Inverses and a certain idea of a "matrix square root" can be computed by applying the desired function to the diagonal. In general, this method of computing the inverse or "square root" is not tremendously efficient (two matrix multiplication operations, the inversion operations, and the spectral decomposition itself) but it can be useful if you already have to compute the spectral decomposition already for some other reason.

*This trick doesn't work as well when matrices are defective or nearly defective: error accumulates to be many times larger than machine precision. This is discussed in detail in the Moler & Van Loan article.

  • 1
    Loosely related, you can also find a "square root" for SPD matrices by essentially the same trick. I was thinking of adding this to my answer, but I think it would fit better here. – Matthew Drury Jul 18 '16 at 15:08
  • @MatthewDrury Yes, that's a good point. I've added some discussion of that. – Sycorax Jul 18 '16 at 15:19
  • @GeneralAbrial Whats the "correct" algorithm for the square root? – Matthew Drury Jul 18 '16 at 15:20
  • 3
    The ability to compute an accurate exponential ultimately tells us that the matrix is diagonalizable (that is, conjugate to a diagonal matrix). That includes far more than the symmetric positive-definite matrices (and far more than the symmetric matrices, for that matter). Thus, although this is a nice property, it wouldn't seem to get to the essence of the matter. – whuber Jul 18 '16 at 16:10
  • 2
    It is a deep theorem of linear algebra that for a symmetric matrix, the Schur decomposition and the spectral decomposition are the same thing. With respect to the matrix square root: there are also iterative methods for it based on Newton-Raphson, but those also work best in the symmetric case. Nevertheless, this answer does not address why SPD matrices are special even among the symmetric matrices. – J. M. is not a statistician Jul 18 '16 at 16:11

There are already several answers explaining why symmetric positive definite matrices are so important, so I will provide an answer explaining why they are not as important as some people, including the authors of some of those answers, think. For the sake of simplicity, I will limit focus to symmetric matrices, and concentrate on Hessians and optimization.

If God had made the world convex, there wouldn't be convex optimization, there would just be optimization. Similarly, there wouldn't be (symmetric) positive definite matrices, there would just be (symmetric) matrices. But that's not the case, so deal with it.

If a Quadratic Programming problem is convex, it can be solved "easily". If it is non-convex, a global optimum can still be found using branch and bound methods (but it may take longer and more memory).

If a Newton method is used for optimization and the Hessian at some iterate is indefinite, then it is not necessary to "finagle" it to positive definiteness. If using a line search, directions of negative curvature can be found and the line search executed along them, and if using a trust region, then there is some small enough trust region such that the solution of the trust region problem achieves descent.

As for Quasi-Newton methods, BFGS (damped if the problem is constrained) and DFP maintain positive definiteness of the Hessian or inverse Hessian approximation. Other Quasi-Newton methods, such as SR1 (Symmetric Rank One) do not necessarily maintain positive definiteness. Before you get all bent out of shape over that, that is a good reason for choosing SR1 for many problems - if the Hessian really isn't positive definite along the path to the optimum, then forcing the Quasi-Newton approximation to be positive definite may result in a lousy quadratic approximation to the objective function. By contrast, the SR1 updating method is "loose as a goose", and can writhely morph its definiteness as it proceeds along.

For nonlinearly constrained optimization problems, what really matters is not the Hessian of the objective function, but the Hessian of the Lagrangian. The Hessian of the Lagrangian may be indefinite even at an (the) optimum, and indeed, it is only the projection of the Hessian of the Lagrangian into the nullspace of the Jacobian of the active (linear and nonlinear) constraints which need be positive semi-definite at the optimum. If you model the Hessian of the Lagrangian via BFGS and thereby constrain it to be positive definite, it might be a terrible fit everywhere, and not work well. By contrast, SR1 can adapt its eigenvalues to what it actually "sees".

There's much more that I could say about all of this, but this is enough to give you a flavor.

Edit: What I wrote 2 paragraphs up is correct. However, I forgot to point out that it also applies to linearly constrained problems. In the case of linearly constrained problems, the Hessian of the Lagrangian is just (reduces down to) the Hessian of the objective function. So the 2nd order optimality condition for a local minimum is that the projection of the Hessian of the objective function into the nullspace of the Jacobian of the active constraints is positive semi-definite. Most notably, the Hessian of the objective function need not (necessarily) be psd at the optimum, and often isn't, even on linearly constrained problems.

You already cited a bunch of reasons why SPD are important yet you still posted the question. So, it seems to me that you need to answer this question first: Why do positive quantities matter?

My answer is that some quantities ought to be positive in order to reconcile with our experiences or models. For instance, the distances between items in the space have to be positive. The coordinates can be negative, but the distances are always non-negative. Hence, if you have a data set and some algorithm that processes it you may well end up with one that breaks down when you feed a negative distance into it. So, you say "my algorithm requires positive distance inputs at all times", and it wouldn't sound like an unreasonable demand.

In the context of statistics, a better analogy would be the variance. So, we calculate the variance as $$\sum_i (x_i-\mu)^2/n$$ It's obvious from the definition that if you feed in the real numbers $x_i$ into the equation the output is always non-negative. Hence, you may build algorithms that work with non-negative numbers, and they may be more efficient than algorithm without this restriction. That's the reason we use them.

So, variance-covariance matrices are positive semi-definite, i.e. "non-negative" in this analogy. The example of an algorithm that requires this condition is Cholesky decomposition, it's very handy. It's often called a "square root of the matrix". So, like the square root of a real number that requires non-negativity, Cholesky wants non-negative matrices. We don't find this constraining when dealing with covariance matrices because they always are.

So, that's my utilitarian answer. The constraints such as non-negativity or SPD allow us build more efficient calculation algorithm or convenient modeling tools that are available when your inputs satisfy these constraints.

Here are two more reasons which haven't been mentioned for why positive-semidefinite matrices are important:

  1. The graph Laplacian matrix is diagonally dominant and thus PSD.

  2. Positive semidefiniteness defines a partial order on the set of symmetric matrices (this is the foundation of semidefinite programming).

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.