Multiplication of large matrix with its transpose

What would be the best way to multiply a large matrix A say, 100000 x 10, by its transpose A' in MATLAB 2021b. The output will be a 100000x100000 matrix. This is giving a memory error > 40 GB.

Previous suggestions in Q&A include avoiding creation of A' before multiplication. This certainly helps in reducing time with smaller matrices as tested with 30,000 points. Thanks.

[Chris Wilson]

Thanks a lot - this is extremely helpful to know - it helped me solve a related problem.

Answer 1

[John D'Errico]

Best is to never create it at all. Instead, use mathematics to work with it. And that depends on what you will do with this matrix. But there are many ways to work with such a matrix, that would completely avoid ever needing to create the matrix at all.

For example, do you want the eigenvalues? I'll do it on a small matrix A, because I want to show how to achieve the same results.

A = randn(1000,5);
[V0,D0] = eig(A*A');
D0(end + [-4:0],end + [-4:0]) % all other eigenvalues will be zero

ans = 5×5

1.0e+03 * 0.9338 0 0 0 0 0 0.9700 0 0 0 0 0 1.0316 0 0 0 0 0 1.0571 0 0 0 0 0 1.0895

% now, without computing A*A' at all.
[Q,R] = qr(A,0);
[V,D] = eig(R*R')

V = 5×5

0.0633 0.2740 -0.3781 -0.7371 0.4845 -0.2697 0.2069 0.8769 -0.3192 0.1168 -0.0450 -0.3803 -0.0440 -0.5893 -0.7100 -0.4761 -0.7565 -0.0318 -0.0119 0.4471 0.8334 -0.4065 0.2919 -0.0859 0.2181

D = 5×5

1.0e+03 * 0.9338 0 0 0 0 0 0.9700 0 0 0 0 0 1.0316 0 0 0 0 0 1.0571 0 0 0 0 0 1.0895

We can also compute the eigenvectors.

My point is, there are many things you can do without ever needing to compute that huge mess of a matrix.

Answer 2

[James Tursa]

The fastest way is to simply write it as A * A', because MATLAB will see that the operands are the same and call a special multi-threaded symmetric BLAS library function to do the multiply without explicitly forming the transpose A' first.

As for your memory problems, that is an entirely different issue. If you really need this result and it is full, then you may be forced to chunk out the operation into smaller pieces and work on them individually. If A has lots of 0's then you could look into sparse representations.