Why are symmetric positive definite (SPD) matrices so important? 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.
 A: 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) = \sqrt{(x-y)^T A (x-y)}
$$
is a metric (also called distance function).   
In addition, positive definite matrices are related to inner 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. 
A: 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.
A: 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.
A: Here are two more reasons which haven't been mentioned for why positive-semidefinite matrices are important:


*

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

*Positive semidefiniteness defines a partial order on the set of symmetric matrices (this is the foundation of semidefinite programming).
A: 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.
A: 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.  
A: 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. 
