Singular value decomposition: Difference between revisions
Content added Content deleted
(Added Wren) |
|||
Line 118: | Line 118: | ||
[[-0.70710678 -0.70710678] |
[[-0.70710678 -0.70710678] |
||
[-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> |
</pre> |