2
$\begingroup$

I am currently using the protocol described in https://stackoverflow.com/questions/55599950/computation-of-pseidoinverse-with-svd-in-c-using-blas-and-lapacke to compute the pseudo inverse of a matrix. Clearly when dealing with large matrices (e.g. 10000 x 10000) the process is extremely expensive and the run becomes quite long. Is it somehow possible to parallelize this process?

Edit: The purpose of the program is to solve a system of equation

$$ A x = y $$

where $A$ is semi-positive definite (i.e. it can be singular). In order to do this I compute the pseudo inverse $A^{-1}$ and compute $x = A^{-1}y$. I am currently using OpenBlas.

$\endgroup$
10
  • 3
    $\begingroup$ 100x100 is very small - not large at all, it should just take a tiny fraction of second and is not worth parallelising. How are you currently trying to do this? If I modify the code in the answer to the question you link to to n = 100 and adjust the set up of a it takes less than 0.1 seconds on my machine. $\endgroup$
    – Ian Bush
    Jan 30, 2023 at 12:37
  • $\begingroup$ Thank you! Sorry for the typo, I was considering 10000x10000. $\endgroup$ Jan 30, 2023 at 14:34
  • 2
    $\begingroup$ Are the "inversions" related in some way? Or are they all independent? Are the matrices dense or sparse? Are they symmetric? Do you have access to a cluster? If not what hardware do you want to use? Could you use a GPU? $\endgroup$
    – Ian Bush
    Jan 30, 2023 at 17:34
  • 2
    $\begingroup$ Do you need an SVD based pseudoinverse for a badly conditioned matrix, or Is the matrix well enough conditioned that you can use the normal equations approach? Do you actually need the pseudoinverse or just the action of the pseudoinverse on one vector (or a small number of vectors)? $\endgroup$ Jan 30, 2023 at 23:51
  • 1
    $\begingroup$ Sorry, that doesn't help very much. Are the A matrices related in any way? Can you simply generate one matrix from the "previous" one? If so how? If not are all all the pseudoinversions independent of one another so we can do may of them at once? $\endgroup$
    – Ian Bush
    Feb 1, 2023 at 12:04

1 Answer 1

3
$\begingroup$

LAPACK, and consequently LAPACKE are technically just interfaces. There are several implementations, some of which are already parallel. OpenBLAS and MKL should scale quite well on multicore machines.You should check which version you link with (I would personally recommend MKL, but OpenBLAS is also fine).

You should also check whether all these operations are really necessary. For example, if the matrix $A$ stays the same each time, you can skip the calculation of the pseudoinverse and only apply it to different vectors.

Calculating the pseudoinverse is also a relatively expensive operation. Do you really need to calculate it? If all you need it for is solving a linear system, then you might save a significant amount of compute time with a QR factorization.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.