# Generating Symmetric Positive Definite Matrices using indices

I was trying to run test cases for CG and I need to generate:

• symmetric positive definite matrices
• of size > 10,000
• FULL DENSE
• Using only matrix indices and if necessary 1 vector (Like $A(i,j) = \dfrac{x(i) - x(j)}{(i+j)}$)

• With condition number less than 1000.

I have tried:

1. Generating random matrices using A=rand(N,N) and then A'A to make it Sym. PD. [This increases condition number]

2. Using the vector appraoch as shown but I can't seem to get a function of (x,i,j) which will ensure Sym and PD.

After much experimentation, I came up with:

a(it,jt) = (vec(it)+vec(jt))/((it-1)^2+(jt-1)^2); If $it \neq jt$

a(it,it) = x(it) if $it=jt$

But this is PD till about 500x500.

1. XLATMR. [With all the grading and scaling, its too difficult to understand. Especially since I cannot understand the underlying linear algebra]

Can someone give me a function in x (vector) and i,j (indices) which meets the above requirements?

• Try adding a large number $\alpha$ (on the order of the condition number) to the diagonal entries of your matrix. This is equivalent to adding $\alpha$ to each of your eigenvalues and should improve the condition number. – Aron Ahmadia Mar 28 '12 at 14:10
• @AronAhmadia, Works brilliantly! Thanks! However, what large number should I add? I tried N itself and it worked till 5000x5000 (just finished simulations), will using a+N*eye(N,N) ensure that it will work for all values beyond 5000? Can you convert your comment into an answer? – Inquest Mar 28 '12 at 17:03

To get a dense positive definite matrix with condition number $c$ cheaply, pick a diagonal matrix $D$ whose diagonal consists of numbers from $[1,c]$ (which will be the eigenvalues), with $1$ and $c$ chosen at least once, and a vector $u$. Then apply a similarity transformation, via Householder transformations, to form the matrix $A:=(I-tuu^T)D(I-tuu^T)$, where $t=2/u^Tu$.

To form this matrix with $O(n^2)$ operations, compute $v:=Du$, $s:=t^2u^Tv/2$, $w:=tv-su$ in $O(n)$ operations and then $A$ as $A=D-uw^T-wu^T$. (If you choose $u$ as the all-one vector, the number of multiplications needed is only $O(n)$.)

Note that the behavior of CG depends a lot on the eigenvalue distribution, which you can easily confirm by varying $D$.

(Adding $\alpha$ to an arbitrary symmetric matrix $A$ needs $\alpha$ larger than the largest eigenvalue. This is hard to compute; so one would have to choose $\alpha=\|A\|$, but then the condition number will generally be very small, not a realistic test case for CG.)

• This is perfect! – Inquest Mar 30 '12 at 16:16
• How should I change the algorithm to provide me matrices when I provide the eigen values? For instance, I want a non-symmetric matrix with eigen values 1,-1,2,-2...50,-50. – Inquest Mar 31 '12 at 9:03
• $D=Diag(-50:50)$ in the above recipe gives you these eigenvalues, but with this choice the matrix is now symmetric indefinite, and CG wouldn't solve such a problem. And in your comment you even ask for a nonsymmetric matrix, wehre CG also doesn't apply. Maybe it is best to pose a new question with what precisely you want. – Arnold Neumaier Mar 31 '12 at 16:42

Try adding a large number $\alpha$ (on the order of the norm of the matrix) to the diagonal entries of your matrix. This is equivalent to adding $\alpha$ to each of your eigenvalues and should improve the condition number by reducing the gap between the largest and smallest eigenvalues.

I'm not sure how you would do it with just one vector, but with two random vectors $x$ and $\theta$ of size $N$, you can produce a positive semi-definite matrix via $$U=\prod_i R_i(\theta_i)\\ A=U\mbox{diag}(\mbox{abs}(x)) U^\ast$$ where $R_i(\cdot)$ is a rotation in the plane of the axes $i$ and $i+1 \mbox{ mod } N$.

If you want to improve the condition number, you can add a fixed positive value to $x$ and rescale if need be.

• I'm sorry but my math isn't all that great yet. Does $U^*$ denote transpose? What does rotation in the plane mean? Also, storing 2 vectors will be a little expensive but still, this look likes a very interesting way. – Inquest Mar 28 '12 at 17:10
• $R_i(\theta_i)=\left(\begin{matrix} 1 & 0 & \cdots \\ 0 & \ddots \\ &&\cos \theta_i & \sin \theta_i \\ & & -\sin \theta_i & \cos \theta_i\\ & & & &1 \\ && & & & \ddots\end{matrix}\right)$ is the rotation matrix. $U^\ast$ denotes the complex conjugate transposed matrix (assuming all is real, it's just the transpose). – Deathbreath Mar 28 '12 at 17:41

An entirely different way to do it would be like this: consider a random vector $x$, then $A=xx^T$ is a rank-one matrix with $N-1$ eigenvalues equal to zero and one strictly positive eigenvalue equal to $\|x\|^2$ with eigenvector $x$. It's also symmetric.

To construct a dense SPD matrix, add together many such rank-1 matrices. In other words, if you have $M$ vectors $x_i$ (e.g. random vectors), then $$A = \sum_{i=1}^M x_i x_i^T$$ is SPD if $M\ge N$ and if the vectors $x_i$ a linearly independent (if they are not linearly independent, then $A$ is positive semidefinite). You can verify if your $M$ vectors (or the first $M\ge N$ vectors you draw from the random process) are linearly independent using Gram-Schmidt successive orthogonalization, but you're also likely to get an SPD matrix if you simply choose $M \gg N$ vectors.