ekmett / ekmett/linear

Singular Value Decomposition

Open
#112 2 comments 1 reaction 0 assignees View on GitHub
feature request
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.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.