Singular Value Decomposition
- Dominant language
- Haskell
- Stars
- 220
- Forks
- 55
- PR merge metrics
- No merged PRs in 30d
Description
#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.)
``` haskell
-- | 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](https://docs.scipy.org/doc/numpy/reference/generated/numpy.diag.html) does the same thing when supplied with a vector for example. Would it be appropriate to open a PR for it?~~
Contributor guide
No contributing guide indexed for this repository
Assessment
This issue has not been assessed yet.