5
$\begingroup$

I'm working with a 5x5 matrix n which has three eigenvalues equal to zero. I want to perform the SVD of this matrix in Mathematica. The matrix n is given by

n[x_] := SetPrecision[{{0, 0, 0, 10.*I, 10.*Sqrt[10./10^x]}, {0, 0, 0, 0.4*10.*I, 0.4*10*Sqrt[10./10^x]}, 
 {0, 0, 0, 0.5*10.*I, 0.5*10*Sqrt[10./10^x]}, {10.*I, 0.4*10.*I, 0.5*10.*I, 10.^x, 0}, 
 {10.*Sqrt[10./10^x], 0.4*10.*Sqrt[10./10^x], 0.5*10.*Sqrt[10./10^x], 0, 10.}}, 1000]; 

The matrix ends up having very small singular values ($10^{-16}$ or so) which I believe are numerical errors. Now consequently, the matrices which come out of SingularValueDecomposition[n[x]] end up being very noisy functions of x.

box1[x_]:=SingularValueDecomposition[n[x],Tolerance->10^-50];

I have tried changing the precision of the input using SetPrecision and the Tolerance of the SVD but nothing appears to remove to noise.

Is there anything else I can do?

Edit: What I do with this after is

box[x_] := SingularValueDecomposition[n[x]]; 
u[x_] := N[Reverse[box[x][[1]], 2] /. x -> 2, 40]; 
\[CapitalLambda][x_] := N[Reverse[box[x][[2]], {1, 2}] /. x -> 2, 40]; 
v[x_] := N[Reverse[box[x][[3]], 2] /. x -> 2, 40]; 
z[x_] := ConjugateTranspose[u[x]] . Conjugate[v[x]]; 
U[x_] := Conjugate[u[x] . MatrixPower[z[x], 1/2]]; 

If I then do DiscretePlot[Abs[u[x][[1, 2]]], {x, -1, 4, 0.01}]:

enter image description here

$\endgroup$
  • $\begingroup$ would you add your Mathematica code please? $\endgroup$ – Mirko Aveta Nov 27 '16 at 14:05
  • $\begingroup$ People here generally like users to post code as Mathematica code instead of images or TeX, so they can copy-paste it. It makes it convenient for them and more likely you will get someone to help you. You may find this this meta Q&A helpful $\endgroup$ – Michael E2 Nov 27 '16 at 14:09
  • $\begingroup$ I've attempted that now although it looks a little unreadable. $\endgroup$ – KMoff Nov 27 '16 at 14:16
  • 1
    $\begingroup$ Don't worry, we can copy and paste it on our own notebooks. $\endgroup$ – Mirko Aveta Nov 27 '16 at 14:17
  • $\begingroup$ Don't you neglect the small singular values (and the associated vectors/matrix columns)? $\endgroup$ – Michael E2 Nov 27 '16 at 14:41
6
$\begingroup$

Is there anything else I can do?

With a matrix of this size it seems better to do the computations symbolically.

Redefine n with exact precision:

n[x_] := {{0, 0, 0, 10*I, 10*Sqrt[10/10^x]}, {0, 0, 0, 4*I, 
    4*Sqrt[10/10^x]}, {0, 0, 0, 5*I, 5*Sqrt[10/10^x]}, {10*I, 4*I, 
    5*I, 10^x, 0}, {10*Sqrt[10/10^x], 4*Sqrt[10/10^x], 
    5*Sqrt[10/10^x], 0, 10}};

n[x] // MatrixForm

enter image description here

Compute SVD symbolically (it is not that slow):

AbsoluteTiming[svdRes = SingularValueDecomposition[n[x]];]

(* {3.53261, Null} *)

Evaluate the result elements at desired precision:

N[svdRes[[1]] /. x -> 2, 40]

(* {{0.09807730004968832306840926757989902604183, 
  0.2925930509160816917625846086760755661284, \
  ...}} *)
$\endgroup$
  • $\begingroup$ Thanks for this. I have updated my question so you can see the noise I get. This remains after performing things symbolically as you suggest. $\endgroup$ – KMoff Nov 27 '16 at 16:26
  • 1
    $\begingroup$ @KMoff It seemed to work for me, if I redefined box1, too: n[x_] := {...<AA's code>...}; box1[x_] = SingularValueDecomposition[n[x]]; -- note the use of = instead of :=. $\endgroup$ – Michael E2 Nov 27 '16 at 17:46
  • $\begingroup$ @MichaelE2 Yes, I have had success with this now in getting rid of the noise from the plots. The issue that is throwing me is that although U[x] now appears to be a smooth function (with x between -2 and 4), the combination of interest to me is Transpose[U[x][[1;;3]].Conjugate[U[x][[1;;3]]]. This still contains noisy entries but I think this may be a separate issue altogether now so perhaps I should accept the answer. $\endgroup$ – KMoff Nov 28 '16 at 12:15
7
$\begingroup$

The problem, as I see it, is that three of the singular values (and eigenvalues) of n[x] are zero, or equivalently, that the 5x5 matrix has rank 2. That means that u[x] and v[x] will carry the first two standard basis vectors to a basis for the column space and orthogonal complement of the null space, respectively. The column vectors of u[x] and v[x] are called the left and right singular vectors, respectively.

When the singular values are distinct and ordered in the traditional decreasing order, the singular vectors are unique up to scalar multiple of absolute value 1. It may be possible to construct the functions u[x] and v[x] to be at least piecewise continuous (by imposing some arbitrary condition, such as some entry in one of the matrices that is generically nonzero be real and nonnegative).

When singular values (or zero eigenvalues) are repeated, then the singular vectors corresponding to a given repeated value need to span a subspace of dimension higher than 1. This makes it harder to consistently choose a basis that varies continuously with the parameter x. When a zero eigenvalue is repeated, as in the OP's case, then the OP is probably right that numerical noise exacerbates the problem. See the last to figures below.

In the OP's example, there are two (nonzero) singular values. These are distinct except for x == 1. The first two columns of u[x] and v[x] vary continuously, except at x == 1. Note however, that the subspaces should be varying continuously with x, I think, although I did not investigate that hypothesis. It seems obvious from the results that SingularValueDecomposition does not yield a continuous parametrization in the OP's form, although @Anton Antonov's symbolic approach seemed to work for me (but not the OP, per OP's comment). Perhaps that will be resolved.

Images, worked out with the original box1[x]. The vector components plotted are flatten out, so that the two singular vectors yield ten components and thus ten plotted curves.

Block[{f},
 f[x_?NumericQ] := Diagonal[box1[x][[2]]];
 ListLinePlot[
  Transpose@Table[Thread@{x, f[x]}, {x, -2, 4, 0.01}], 
  PlotLabel -> "Singular values"]
 ]

Mathematica graphics

Block[{f},
 f[x_?NumericQ] := Flatten[box1[x][[1, All, {1, 2}]]];
 Module[{n = 1}, 
  Plot[f[x], {x, -2, 4}, 
    PlotLabel -> "First two left singular vector components"] /. 
   l_Line :> {ColorData[97, n++], l}]
 ]

Mathematica graphics

Block[{f},
 f[x_?NumericQ] := Flatten[box1[x][[3, All, {1, 2}]]];
 Module[{n = 1}, 
  Plot[f[x], {x, -2, 4}, 
    PlotLabel -> "First two right singular vector components"] /. 
   l_Line :> {ColorData[97, n++], l}]
 ]

Mathematica graphics

Block[{f},
 f[x_?NumericQ] := Flatten[box1[x][[3, All, {3, 4, 5}]]];
 Module[{n = 1}, 
  Plot[Re@f[x], {x, -2, 4}, 
    PlotLabel -> "Last three right singular vector components"] /. 
   l_Line :> {ColorData[97, n++], l}]
 ]

Mathematica graphics

DiscretePlot[Im@Det@box1[x][[1]], {x, -1, 4, 0.01}, Frame -> True, 
 PlotLabel -> "Imaginary part of Det[u[x]] varies wildly"]

Mathematica graphics

$\endgroup$
  • $\begingroup$ Thanks very much for this very informative answer. I have one question - when you say "It may be possible to construct the functions u[x] and v[x] to be at least piecewise continuous (by imposing some arbitrary condition, such as some entry in one of the matrices that is generically nonzero be real and nonnegative)." what did you have in mind exactly? $\endgroup$ – KMoff Nov 30 '16 at 18:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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