Generate two matrices such that multiplication of these matrices results a symmetric positive definite matrix

Question

I want to generate a tall&skinny matrix and a fat&short matrix, such that multiplication of these matrices results a symmetric positive definite matrix.

My inefficient/expensive solution is:

A       = rand(5,5);
A       = A + A';
P       = A + n*eye(n);
% P is a symmetric positive definite matrix 
[U,S,V] = svd(P);
X       = U * sqrt(S);
X2      = sqrt(S) * V';

X and X2 contain the desired matrices. I do not want to run an expensive operation such as svd(). I want to generate two matrices as cheap as possible, i.e., the solution must be fast and must avoid constructing the symmetric positive definite matrix.

The first matrix (X) must be tall&skinny. The second matrix (X2) must be fat&short.

Answer 1

You're right, that's pretty expensive. You cannot make a positive definite matrix by just multiplying any skinny * fat. Instead, you can do one of the following: skinny' * skinny or fat * fat' + c*identity where c > 0 where skinny and fat in the above are fixed (so you're multiply a matrix by its own transpose.

Here's a faster way:

m = 5; n = 3; % these can be any natural number as long as n<m
A = rand(m,n); % a tall skinny matrix

posDef = A'*A % will be positive definite n*n matrix with probability 1

or you can replace that last line with

lambdaMin = 1 % any positive number that you want to be your smallest eigenvalue
posDef = A*A' + lambdaMin*eye(m) % will be positive definite m*m matrix with probability 1

Why this works. Any finite-dimensional, real positive definite matrix can be decomposed into a sum of dyads of the form $a \otimes a = a a^\top$ (an outer product of a finite-dimensional, real vector and itself). These dyads have rank one, are symmetric, and are in fact positive semidefinite. What we're doing in the code above is just stacking these vectors side by side in a matrix $A$ so that $A^\top A = \sum_{i=1}^m a_i a_i^\top$ where $a_i$ is the $i^{th}$ row of $A$. So long as no one dyads is in the span of the other dyads (which will happen almost surely when we randomly generate $A$), the sum of the dyads remains symmetric and has either full-rank or rank as large as the number of elements of the sum (whichever rank is smaller).