# Generate symmetric positive definite matrix with a pre-specified sparsity pattern

I am trying to generate a correlation matrix $p\times p$ (symmetric p.s.d) with a pre-specified sparsity structure (specified by a graph on $p$ nodes). The nodes that are connected in the graph have correlation $\rho \sim U(0,1)$, rest all are 0 and the diagonal is all 1.

I have tried generating this matrix several times but only rarely get a valid correlation matrix.

Is there a way that I can assure a correlation matrix w.h.p? Note that I can only have positive correlation so $\rho \sim U(-1,1)$ etc. is not an option.

Any help is greatly appreciated!

• Maybe the function nearPD of the package Matrix in R can help. – niandra82 Jan 1 '15 at 12:53
• What is your measure of sparsity that is fixed for you? Should your data be binary or nonnegative continuous? – ttnphns Jan 1 '15 at 13:07
• @niandra82: nearPD is no good as it will destroy the sparsity of the matrix. – Blade Runner Jan 1 '15 at 14:55
• In general there are no such matrix distributions as described in this question. Consider, for instance, the $3\times 3$ case with three coefficients $\rho,\sigma,\tau$. If $\tau=0$ and $\rho\gt 0, \sigma\gt 0$, then $\rho^2+\sigma^2\lt 1$ if and only if the matrix is positive definite. But then you cannot have both $\rho\sim U(0,1)$ and $\sigma\sim U(0,1)$. – whuber Jan 1 '15 at 20:20
• Then why not generate the correlation matrix first. Then create a symetric index for that matrix where you force the indexed elements to 0. Sparcity would be specified by the size of the index and you can incorperate randommess through a function like sample in r. No matter how many off diagonal elements you force to 0, the matix will still be p.d. – Zachary Blumenfeld Jan 1 '15 at 21:03

Close but no cigar for @Rodrigo de Azevedo.

The solution is to use semidefinite programming to find the maximum value, $\rho_{max}$, and the minimum value (subject to being nonnegative), $\rho_{min}$, of $\rho$ such that the correlation matrix with prescribed sparsity pattern is positive semidefinite (psd). All values of $\rho$ such that $\rho_{max}\le \rho \le \rho_{max}$, will produce psd matrices (exercise for reader)

Therefore, you must either choose a distribution of $\rho$ which can only take values in $[\rho_{max},\rho_{max}]$, or you must use acceptance/rejection and reject any generated values of $\rho$ which do not produce a psd matrix.

Example for a 4 by 4 matrix using YALMIP under MATLAB

sdpvar rho % declare rho to be a scalar variable
% find maximum value of rho (by minimizing -rho) subject to prescribed matrix being psd.
optimize([1 0 rho 0;0 1 rho 0;rho rho 1 rho;0 0 rho 1] >= 0,-rho)
% find minimum value of rho subject to prescribed matrix being psd and rho being >= 0.
optimize([[1 0 rho 0;0 1 rho 0;rho rho 1 rho;0 0 rho 1] >= 0,rho >= 0],rho)


Results: maximum rho = 0.57735, minimum rho = 0. It is readily apparent that zero will be the minimum value of rho subject to rho being nonegative and the prescribed matrix being psd, regardless of dimension or sparsity pattern. Therefore, it is unnecessary to run the semidefinite optimization to find the minimum nonnegative value of $\rho$.

• This is an interesting interpretation of the question: it assumes all nonzero off-diagonal coefficients are equal (thereby tremendously simplifying the problem). It's not clear whether that was the intended interpretation, or whether all nonzero off-diagonal coefficients are supposed to be independent realizations from a common distribution. – whuber Jan 30 '17 at 22:22
• That is the interpretation I made. Now that you mention it, I could see a different interpretation being possible. At least my interpretation has the virtue of resulting in a fairly well-defined problem. I suppose a problem can be formulated, whose solution I have not researched, to find the maximum value of ρ such that all non-zero off-diagonal elements of one triangle of the correlation matrix can be filled in with not necessarily equal nonnegative values ≤ that value, and necessarily make the fully populated matrix be psd. – Mark L. Stone Jan 30 '17 at 22:37

A correlation matrix is symmetric, positive semidefinite, and has $1$'s on its main diagonal. One can find a $n \times n$ correlation matrix by solving the following semidefinite program (SDP) where the objective function is arbitrary, say, the zero function

$$\begin{array}{ll} \text{minimize} & \mathrm \langle \mathrm O_n, \mathrm X \rangle\\ \text{subject to} & x_{11} = x_{22} = \cdots = x_{nn} = 1\\ & \mathrm X \succeq \mathrm O_n\end{array}$$

If one has additional constraints, such as sparsity constraints

$$x_{ij} = 0 \text{ for all } (i,j) \in \mathcal Z \subset [n] \times [n]$$

and non-negativity constraints, $\mathrm X \geq \mathrm O_n$, then one solves the following SDP

$$\begin{array}{ll} \text{minimize} & \mathrm \langle \mathrm O_n, \mathrm X \rangle\\ \text{subject to} & x_{11} = x_{22} = \cdots = x_{nn} = 1\\ & x_{ij} = 0 \text{ for all } (i,j) \in \mathcal Z \subset [n] \times [n]\\ & \mathrm X \geq \mathrm O_n\\ & \mathrm X \succeq \mathrm O_n\end{array}$$

A $3 \times 3$ example

Suppose we want to have $x_{13} = 0$ and $x_{12}, x_{23} \geq 0$. Here's a MATLAB + CVX script,

cvx_begin sdp

variable X(3,3) symmetric

minimize( trace(zeros(3,3)*X) )
subject to

% put ones on the main diagonal
X(1,1)==1
X(2,2)==1
X(3,3)==1

% put a zero in the northeast and southwest corners
X(1,3)==0

% impose nonnegativity
X(1,2)>=0
X(2,3)>=0

% impose positive semidefiniteness
X >= 0

cvx_end


Running the script,

Calling sedumi: 8 variables, 6 equality constraints
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 6, order n = 6, dim = 12, blocks = 2
nnz(A) = 8 + 0, nnz(ADA) = 36, nnz(L) = 21
it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
0 :            3.00E+000 0.000
1 : -1.18E-001 6.45E-001 0.000 0.2150 0.9000 0.9000   1.86  1  1  1.2E+000
2 : -6.89E-004 2.25E-002 0.000 0.0349 0.9900 0.9900   1.52  1  1  3.5E-001
3 : -6.48E-009 9.72E-007 0.097 0.0000 1.0000 1.0000   1.01  1  1  3.8E-006
4 : -3.05E-010 2.15E-009 0.000 0.0022 0.9990 0.9990   1.00  1  1  1.5E-007
5 : -2.93E-016 5.06E-015 0.000 0.0000 1.0000 1.0000   1.00  1  1  3.2E-013

iter seconds digits       c*x               b*y
5      0.3   5.8  0.0000000000e+000 -2.9302886987e-016
|Ax-b| =  1.7e-015, [Ay-c]_+ =  6.1E-016, |x|= 2.0e+000, |y|= 1.5e-015

Detailed timing (sec)
Pre          IPM          Post
1.563E-001    2.500E-001    1.094E-001
Max-norms: ||b||=1, ||c|| = 0,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 1.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +0


Let's see what solution CVX found,

>> X

X =

1.0000    0.4143         0
0.4143    1.0000    0.4143
0    0.4143    1.0000


Is this matrix positive semidefinite? Positive definite?

>> rank(X)

ans =

3

>> eigs(X)

ans =

1.5860
1.0000
0.4140


It is positive definite, as expected. We can find positive semidefinite correlation matrices by choosing a nonzero (linear) objective function.

• Because on this site "generate" would be understood to mean "draw from a random distribution," could you explain how your code produces random correlation matrices and indicate what distribution they follow? – whuber Oct 18 '16 at 13:47
• @whuber The OP is asking for the impossible. You commented on that on January 1, 2015. If you want to generate random correlation matrices, then generate a random square matrix and use it in the objective function in the semidefinite program above. Or, generate realizations of a random variable that is uniform over the cube $$[-1,1]^{\binom{n}{2}}$$ put them in the off-diagonal entries of (correlation) matrices with $1$'s on the main diagonal, and discard the ones that aren't positive semidefinite. If there are nonnegativity constraints, then uniformly sample the cube $$[0,1]^{\binom{n}{2}}$$ – Rodrigo de Azevedo Oct 18 '16 at 13:55
• @whuber Here's the 3D elliptope [png], which maps to the set of $3 \times 3$ correlation matrices. What the OP wants is to intersect the elliptope with the nonnegative octant, then intersect it with planes of the form $x_{ij} = 0$. If the matrix is $\succ 0$, then it must be in the interior of the elliptope. Using SDP with nonzero objective functions, one can sample the surface of the elliptope. As the elliptope is convex, convex combinations of surface points will also map to correlation matrices. – Rodrigo de Azevedo Oct 18 '16 at 14:14
• That's an excellent way to describe the situation. – whuber Oct 18 '16 at 14:20
• You are correct about how the relative volumes shrink. That is precisely why this is a difficult problem. – whuber Oct 18 '16 at 16:09