4
$\begingroup$

The following is a small part of my research:

I got a 4 by 4 matrix $\mathbf{M}$ of rank 3:

M={{1, 1, 1, 1}, {2, -2, 2.0006219879401477, -2.0006219879401477}, {4.193503236170054, 0.07680669095242888, 4.195764689057942, 0.0767525731461281}, {19.781537502364305, -0.36231149755725733, 19.793967330980973, -0.36208844823571257}}

I need to solve the system $\mathbf{MC=0}$ where $\mathbf{C}=[C_1,C_2,C_3,C_4]^\top$. I know that $\mathbf{M}$ has rank equal to 3, but due to computation error, Mathematica gives its rank equal to 4. (One can check using the command "Eigenvalues[M]" to find that one of its eigenvalues is very close to zero.)

I set $C_1=1$, and want to solve for $C_2,C_3,C_4$. I'm wondering if there is a way I can get a good approximation of $C_2,C_3,C_4$. I thought of using Gaussian Elimination, but I have no idea how to proceed, as I cannot change the fact that the rank is 4 in Mathemtaica.

Any help or hint will be much appreciated!

$\endgroup$
3
  • 1
    $\begingroup$ The Tolerance option ofNullSpace is good for this situation. In[2220]:= mat = {{1, 1, 1, 1}, {2, -2, 2.0006219879401477, -2.0006219879401477}, {4.193503236170054, 0.07680669095242888, 4.195764689057942, 0.0767525731461281}, {19.781537502364305, -0.36231149755725733, 19.793967330980973, -0.36208844823571257}}; NullSpace[mat, Tolerance -> 10^(-6)] Out[2221]= {{0.23992, -0.665256, -0.23978, 0.665115}} $\endgroup$ Dec 20, 2023 at 22:02
  • $\begingroup$ @DanielLichtblau NullSpace gives exactly the same result as a SingularValueDecomposition (see update in my solution) and seems much simpler! $\endgroup$
    – Roman
    Dec 21, 2023 at 7:32
  • $\begingroup$ @Roman As you may have guessed, indeed NullSpace for approximate values does use the SVD. $\endgroup$ Dec 21, 2023 at 14:57

1 Answer 1

6
$\begingroup$
M = {{1, 1, 1, 1},
     {2, -2, 2.0006219879401477, -2.0006219879401477},
     {4.193503236170054, 0.07680669095242888, 4.195764689057942, 0.0767525731461281},
     {19.781537502364305, -0.36231149755725733, 19.793967330980973, -0.36208844823571257}};

Look for the eigenvector associated with the smallest eigenvalue (by absolute value):

c = MinimalBy[Transpose[Eigensystem[M]], Abs@*First][[1, 2]]
(*    {0.240408, -0.66508, -0.240268, 0.664939}    *)

The result is pretty close to $M\cdot c=0$:

M . c
(*    {-3.98559*10^-7, 1.1026*10^-6, 3.98325*10^-7, -1.10236*10^-6}    *)

Normalizing by the first element $c_1=1$:

cn = Normalize[c, First]
(*    {1., -2.76646, -0.999414, 2.76587}    *)

If your matrix is large, you could speed things up by trying to compute only the smallest eigenvalue/eigenvector pair:

Eigensystem[M, -1]
(*    {{-1.65784*10^-6}, {{0.240408, -0.66508, -0.240268, 0.664939}}}    *)

and therefore

c = Eigenvectors[M, -1][[1]]
(*    {0.240408, -0.66508, -0.240268, 0.664939}    *)

update: singular-value decomposition

On second thought, an eigenvalue decomposition may not be the best way forward and a more general singular-value decomposition could give a slightly more accurate result.

The singular value decomposition is

{u, σ, v} = SingularValueDecomposition[M];

Look for the smallest singular value by absolute magnitude:

s = OrderingBy[Diagonal[σ], Abs, 1][[1]]
(*    4    *)

The minimizing vector is the corresponding vector in $v$:

d = Transpose[v][[s]]
(*    {0.23992, -0.665256, -0.23978, 0.665115}    *)

See that it fits better:

{Norm[M . c], Norm[M . d]}
(*    {1.65784*10^-6, 5.37312*10^-7}    *)

Faster by computing only the smallest singular value by Arnoldi iteration:

d = SingularValueDecomposition[M, -1,
      Method -> {"Arnoldi",
                 "Criteria" -> "Magnitude",
                 "MaxIterations" -> 10^6}][[3, All, 1]]
(*    {-0.23992, 0.665256, 0.23978, -0.665115}    *)

Even faster and fully equivalent: compute the smallest eigenvalue and its associated eigenvector of $M^T\cdot M$: this is the vector that minimizes $\|M\cdot\vec{d}\|^2=\vec{d}\cdot (M^T\cdot M)\cdot \vec{d}$ under the constraint that $\|\vec{d}\|^2=\vec{d}\cdot\vec{d}=1$:

d = Eigenvectors[Transpose[M] . M, -1,
      Method -> {"Arnoldi",
                 "Criteria" -> "Magnitude",
                 "MaxIterations" -> 10^6}][[1]]
(*    {-0.23992, 0.665256, 0.23978, -0.665115}    *)
$\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.