Singular value decomposition: Difference between revisions
Added Wren
(Added Wren) |
|||
Line 118:
[[-0.70710678 -0.70710678]
[-0.70710678 0.70710678]]
</pre>
=={{header|Wren}}==
{{libheader|Wren-matrix}}
{{libheader|Wren-fmt}}
We don't have a built-in method for this so we need to work from first principles.
The following only decomposes 2 x 2 matrices.
<syntaxhighlight lang="ecmascript">import "./matrix" for Matrix
import "./fmt" for Fmt
var svd2x2 = Fn.new { |a, width, prec|
var at = a.transpose
var p = a * at
// the eigenvalues, λ, are the solutions of the characteristic polynomial:
// λ^2 - (p[0,0] + p[1,1]) * λ + p[0,0]*p[1,1] - p[0,1]*p[1,0] = 0
var b = -(p[0, 0] + p[1, 1])
var c = p[0, 0]*p[1, 1] - p[0, 1]*p[1, 0]
var d = (b * b - 4 * c).sqrt
var l1 = (-b + d) / 2
var l2 = (-b - d) / 2
// hence the singular values and 'sigma' are:
var s1 = l1.sqrt
var s2 = l2.sqrt
var sigma = Matrix.new([[s1, 0], [0, s2]])
// now find the eigenvectors
p = at * a
var i = Matrix.identity(p.numRows)
var m1 = p - i*l1
var m2 = p - i*l2
m1.toReducedRowEchelonForm
m2.toReducedRowEchelonForm
// solve for null space and convert to unit vector
var ev1 = Matrix.new([[-m1[0,0] / m1[0, 1]], [1]])
ev1 = ev1 / ev1.norm
var ev2 = Matrix.new([[-m2[0,0] / m2[0, 1]], [1]])
ev2 = ev2 / ev2.norm
// we can now put these vectors together to get 'v' transpose
var vt = Matrix.new([[ev1[0, 0], ev2[0, 0]], [ev1[1, 0], ev2[1, 0]]])
// and since 'v' is orthogonal its transpose and inverse are the same
var u = a * vt * (sigma.inverse)
System.print("U:")
Fmt.mprint(u, width, prec)
System.print("\nΣ:")
Fmt.mprint(sigma, width, prec)
System.print("\nVT:")
Fmt.mprint(vt, width, prec)
}
var a = Matrix.new([[3, 0], [4, 5]])
svd2x2.call(a, 17, 14)</syntaxhighlight>
{{out}}
<pre>
U:
| 0.31622776601684 -0.94868329805051|
| 0.94868329805051 0.31622776601684|
Σ:
| 6.70820393249937 0.00000000000000|
| 0.00000000000000 2.23606797749979|
VT:
| 0.70710678118655 -0.70710678118655|
| 0.70710678118655 0.70710678118655|
</pre>
| |||