3
$\begingroup$

I am working with Conjugate Gradient method to solve for $Ax = b$, where $A$ is an extremely large PSD and Singular matrix. I cannot directly access the elements of $A$, but I can compute $Av$ for any vector $v$.

Now, I would like to solve for $Ax = 0$, which is effectively finding one vector in the null space of $A$. I understand that it is slow to solve it with CG directly. I wonder if there are other iterative methods that would give a reasonably good solution with only the access of $Av$ very quickly?

$\endgroup$
6
  • $\begingroup$ Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, $\endgroup$
    – Dunlop
    Mar 2, 2023 at 4:43
  • $\begingroup$ 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. $\endgroup$
    – Dunlop
    Mar 2, 2023 at 4:43
  • 1
    $\begingroup$ does this question have anything to do with the software Mathematica? If not it might be better at other Stack Exchange sites, and may risk closure $\endgroup$
    – Dunlop
    Mar 2, 2023 at 4:44
  • $\begingroup$ @Dunlop this question has a clear mathematical answer: use the Arnoldi iteration, as implemented in ARPACK, to find the eigenvalue/vector with smallest magnitude. There is no mathematical difficulty here, it's all solved and well-known. Unfortunately, ARPACK's interface that only requires the black-box operation $\vec{v}\mapsto A\cdot\vec{v}$ is not exposed in Mathematica, and so the question becomes a particular difficulty in Mathematica. $\endgroup$
    – Roman
    Mar 2, 2023 at 11:24
  • $\begingroup$ @Roman that makes things clearer. It was just the question didnt seem as though it had anything to do with Mathematica $\endgroup$
    – Dunlop
    Mar 2, 2023 at 11:26

2 Answers 2

4
$\begingroup$

I'm implementing a simplified version of the Arnoldi iteration with spectral shift here.

First, create a random 10✕10 PSD matrix that we'll use as an example:

n = 10;
A = Transpose[#] . # &@RandomVariate[NormalDistribution[], {n, n}];

The eigenvalues of this matrix are all nonnegative (in your case, some will even be zero):

Eigenvalues[A] // Sort
(*    {0.0540622, 0.321272, 0.680714, 2.80777, 4.01219,
       10.5359, 12.1004, 17.871, 21.7634, 34.2863}         *)

Find the largest eigenvalue/vector by magnitude: a straight-forward Arnoldi iteration,

vmax = FixedPoint[Normalize[A . #] &, RandomReal[{0, 1}, n],
                  SameTest -> (Norm[#1 - #2] < 10^-12 &)]
(*    {-0.334135, 0.233949, 0.424574, -0.437277, 0.115131,
       0.376451, -0.0107379, 0.471395, 0.0544037, 0.286161}    *)
emax = (A . vmax) . vmax
(*    34.2863    *)

Now we use a "spectral shift" operation: we set the matrix $B=\lambda_{\text{max}}I-A$, whose largest eigenvalue is $\lambda_{\text{max}}-\lambda_{\text{min}}$. Therefore, finding the largest eigenvalue of $B$ allows us to find the smallest eigenvalue/vector of $A$ by magnitude:

vmin = FixedPoint[Normalize[emax # - A . #] &, RandomReal[{0, 1}, n],
                  SameTest -> (Norm[#1 - #2] < 10^-12 &)]
(*    {0.249744, -0.562473, 0.237567, -0.211816, 0.294718,
       -0.246993, -0.0820622, 0.0486029, -0.52359, 0.298061}    *)
emin = (A . vmin) . vmin
(*    0.0540622    *)

We see that emin and emax are the correct values (as calculated with Eigenvalues[A] above). We've found them iteratively, only requiring the operation $\vec{v}\mapsto A\cdot\vec{v}$ without ever looking at $A$ explicitly. You can replace every mention of A . # in the above code with a black-box function Aop[#]:

Aop[v_ /; VectorQ[v, NumericQ]] := A . v

vmax = FixedPoint[Normalize[Aop[#]] &, RandomReal[{0, 1}, n],
                  SameTest -> (Norm[#1 - #2] < 10^-12 &)];
emax = Aop[vmax] . vmax;
vmin = FixedPoint[Normalize[emax # - Aop[#]] &, RandomReal[{0, 1}, n],
                  SameTest -> (Norm[#1 - #2] < 10^-12 &)]
emin = Aop[vmin] . vmin

(*    {0.249744, -0.562473, 0.237567, -0.211816, 0.294718,
       -0.246993, -0.0820622, 0.0486029, -0.52359, 0.298061}    *)

(*    0.0540622    *)

If your matrix is singular, then the smallest eigenvalue emin will be zero; but this assumption is not made in the above code.

$\endgroup$
8
  • $\begingroup$ Thank you very much for your answer. However, your solution seems to require access to eigenvalues of A, which is also inaccessible. $\endgroup$ Mar 2, 2023 at 12:20
  • $\begingroup$ No it doesn't. I only calculated the eigenvalues to show that the iterative results are correct. $\endgroup$
    – Roman
    Mar 2, 2023 at 12:55
  • $\begingroup$ Oh I see. Essentially you compute emax with Arnoldi iteration, then use emax to find vmin. This is indeed a method that fits the question description. However, this method seems a little bit time-consuming when Av operation is expensive, correct? $\endgroup$ Mar 2, 2023 at 17:49
  • $\begingroup$ It depends on the SameTest precision that you need. If a sloppy result is good enough, this method can be very fast. Also, if you have an a priori estimate of an upper bound of emax, then you can use that directly without calculating it by iteration. $\endgroup$
    – Roman
    Mar 2, 2023 at 18:54
  • $\begingroup$ Do you have any recommendation on how to acquire an accurate priori estimate? $\endgroup$ Mar 2, 2023 at 21:08
3
$\begingroup$

Provided that A x can be calculated fast, you could try to minimize the norm of Ax where x is unit vector. E.g.:

For a test we create a singular matrix by:

n = 20;
m = DiagonalMatrix[Append[RandomReal[{-1, 1}, n - 1], 0]];
mo = Orthogonalize[RandomReal[{-1, 1}, {n, n}]];
m = Transpose[mo] . m . mo;
Eigenvalues[m]

enter image description here

We now minimize Ax with Norm[x]==1:

vs = Array[Subscript[v, #] &, n]; sol = vs /. Minimize[{Norm[m . vs], Norm[vs] == 1}, vs][[2]]

To test if sol is in the Nullspace:

m . sol // Abs // Max
3.44922*10^-7
$\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.