Singular Value Decomposition

[rehno-lindeque]

86 already mentions adding SVD, but I thought it deserved a separate issue to track progress and would be very convenient for calculating eigenvectors / eigenvalues.

I have a largely untested and unbenchmarked implementation of SVD for 2x2 matrices I'm working on that might be of some help to someone who doesn't feel like pulling in more dependencies and marshaling between them. It's likely that this one isn't the most optimal implementation, I've lifted it from http://scicomp.stackexchange.com/a/19646/3470. (I'd appreciate references to other 2x2 SVD implementations if anyone has suggestions.)

-- | Singular value decomposition on a 2x2 matrix
svd22 :: M22 Double -> (M22 Double, V2 Double, M22 Double)
svd22 m@(V2 (V2 m00 m01) (V2 m10 m11)) =
  let v1@(V2 x1 y1) = V2 (m00 - m11) (m10 + m01)
      v2@(V2 x2 y2) = V2 (m00 + m11) (m10 - m01)

      h1 = (sqrt . sum . sqr) v1 -- hypotenuse (vector length)
      h2 = (sqrt . sum . sqr) v2

      t1 = x1 / h1
      t2 = x2 / h2

      cc = sqrt ((1 + t1) * (1 + t2))
      ss = sqrt ((1 - t1) * (1 - t2))
      cs = sqrt ((1 + t1) * (1 - t2))
      sc = sqrt ((1 - t1) * (1 + t2))

      (c1,s1) = ((cc - ss) / 2, (sc + cs) / 2)
      u1 = V2 (V2 c1 (-s1)) (V2 s1 c1)

      sigma@(V2 sigma1 sigma2) = V2 ((h1 + h2) / 2) ((h1 - h2) / 2)

      u2 = scaled (V2 (1 / sigma1) (if h1 /= h2 then 1 / sigma2 else 0)) !*! transpose u1 !*! m
  in (u1, sigma, u2)
  where
    sqr v =  v * v

As an aside, I believe fromDiagonal might be a helpful utility function. Numpy's diag does the same thing when supplied with a vector for example. Would it be appropriate to open a PR for it?

[acowley]

You can use scaled for fromDiagonal. I haven't yet looked closely at this code, but I'm in support of the idea.

[rehno-lindeque]

Thank you, I updated my example code with scaled.