2f:T1ecf,/** SVD procedure as explained in "Singular Value Decomposition and Least Squares Solutions. By G.H. Golub et al." * * This procedure computes the singular values and complete orthogonal decomposition of a real rectangular matrix A: * A = U * diag(q) * V(t), U(t) * U = V(t) * V = I * where the arrays a, u, v, q represent A, U, V, q respectively. The actual parameters corresponding to a, u, v may * all be identical unless withu = withv = {true}. In this case, the actual parameters corresponding to u and v must * differ. m >= n is assumed (with m = a.length and n = a[0].length) * * @param a {Array} Represents the matrix A to be decomposed * @param [withu] {bool} {true} if U is desired {false} otherwise * @param [withv] {bool} {true} if U is desired {false} otherwise * @param [eps] {Number} A constant used in the test for convergence; should not be smaller than the machine precision * @param [tol] {Number} A machine dependent constant which should be set equal to B/eps0 where B is the smallest * positive number representable in the computer * * @returns {Object} An object containing: * q: A vector holding the singular values of A; they are non-negative but not necessarily ordered in * decreasing sequence * u: Represents the matrix U with orthonormalized columns (if withu is {true} otherwise u is used as * a working storage) * v: Represents the orthogonal matrix V (if withv is {true}, otherwise v is not used) * */ SVD3 = (a, withu, withv, eps, tol) => { // Define default parameters withu = withu !== undefined ? withu : true withv = withv !== undefined ? withv : true eps = eps || Math.pow(2, -52) tol = 1e-64 / eps // throw error if a is not defined if (!a) { throw new TypeError('Matrix a is not defined') } // Householder's reduction to bidiagonal form const n = a[0].length const m = a.length if (m < n) { throw new TypeError('Invalid matrix: m < n') } let i, j, k, l, l1, c, f, g, h, s, x, y, z g = 0 x = 0 const e = [] const u = [] const v = [] const mOrN = (withu === 'f') ? m : n // Initialize u for (i = 0; i < m; i++) { u[i] = new Array(mOrN).fill(0) } // Initialize v for (i = 0; i < n; i++) { v[i] = new Array(n).fill(0) } // Initialize q const q = new Array(n).fill(0) // Copy array a in u for (i = 0; i < m; i++) { for (j = 0; j < n; j++) { u[i][j] = a[i][j] } } for (i = 0; i < n; i++) { e[i] = g s = 0 l = i + 1 for (j = i; j < m; j++) { s += Math.pow(u[j][i], 2) } if (s < tol) { g = 0 } else { f = u[i][i] g = f < 0 ? Math.sqrt(s) : -Math.sqrt(s) h = f * g - s u[i][i] = f - g for (j = l; j < n; j++) { s = 0 for (k = i; k < m; k++) { s += u[k][i] * u[k][j] } f = s / h for (k = i; k < m; k++) { u[k][j] = u[k][j] + f * u[k][i] } } } q[i] = g s = 0 for (j = l; j < n; j++) { s += Math.pow(u[i][j], 2) } if (s < tol) { g = 0 } else { f = u[i][i + 1] g = f < 0 ? Math.sqrt(s) : -Math.sqrt(s) h = f * g - s u[i][i + 1] = f - g for (j = l; j < n; j++) { e[j] = u[i][j] / h } for (j = l; j < m; j++) { s = 0 for (k = l; k < n; k++) { s += u[j][k] * u[i][k] } for (k = l; k < n; k++) { u[j][k] = u[j][k] + s * e[k] } } } y = Math.abs(q[i]) + Math.abs(e[i]) if (y > x) { x = y } } // Accumulation of right-hand transformations if (withv) { for (i = n - 1; i >= 0; i--) { if (g !== 0) { h = u[i][i + 1] * g for (j = l; j < n; j++) { v[j][i] = u[i][j] / h } for (j = l; j < n; j++) { s = 0 for (k = l; k < n; k++) { s += u[i][k] * v[k][j] } for (k = l; k < n; k++) { v[k][j] = v[k][j] + s * v[k][i] } } } for (j = l; j < n; j++) { v[i][j] = 0 v[j][i] = 0 } v[i][i] = 1 g = e[i] l = i } } // Accumulation of left-hand transformations if (withu) { if (withu === 'f') { for (i = n; i < m; i++) { for (j = n; j < m; j++) { u[i][j] = 0 } u[i][i] = 1 } } for (i = n - 1; i >= 0; i--) { l = i + 1 g = q[i] for (j = l; j < mOrN; j++) { u[i][j] = 0 } if (g !== 0) { h = u[i][i] * g for (j = l; j < mOrN; j++) { s = 0 for (k = l; k < m; k++) { s += u[k][i] * u[k][j] } f = s / h for (k = i; k < m; k++) { u[k][j] = u[k][j] + f * u[k][i] } } for (j = i; j < m; j++) { u[j][i] = u[j][i] / g } } else { for (j = i; j < m; j++) { u[j][i] = 0 } } u[i][i] = u[i][i] + 1 } } // Diagonalization of the bidiagonal form eps = eps * x let testConvergence for (k = n - 1; k >= 0; k--) { for (let iteration = 0; iteration < 50; iteration++) { // test-f-splitting testConvergence = false for (l = k; l >= 0; l--) { if (Math.abs(e[l]) <= eps) { testConvergence = true break } if (Math.abs(q[l - 1]) <= eps) { break } } if (!testConvergence) { // cancellation of e[l] if l>0 c = 0 s = 1 l1 = l - 1 for (i = l; i < k + 1; i++) { f = s * e[i] e[i] = c * e[i] if (Math.abs(f) <= eps) { break // goto test-f-convergence } g = q[i] q[i] = Math.sqrt(f * f + g * g) h = q[i] c = g / h s = -f / h if (withu) { for (j = 0; j < m; j++) { y = u[j][l1] z = u[j][i] u[j][l1] = y * c + (z * s) u[j][i] = -y * s + (z * c) } } } } // test f convergence z = q[k] if (l === k) { // convergence if (z < 0) { // q[k] is made non-negative q[k] = -z if (withv) { for (j = 0; j < n; j++) { v[j][k] = -v[j][k] } } } break // break out of iteration loop and move on to next k value } // Shift from bottom 2x2 minor x = q[l] y = q[k - 1] g = e[k - 1] h = e[k] f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2 * h * y) g = Math.sqrt(f * f + 1) f = ((x - z) * (x + z) + h * (y / (f < 0 ? (f - g) : (f + g)) - h)) / x // Next QR transformation c = 1 s = 1 for (i = l + 1; i < k + 1; i++) { g = e[i] y = q[i] h = s * g g = c * g z = Math.sqrt(f * f + h * h) e[i - 1] = z c = f / z s = h / z f = x * c + g * s g = -x * s + g * c h = y * s y = y * c if (withv) { for (j = 0; j < n; j++) { x = v[j][i - 1] z = v[j][i] v[j][i - 1] = x * c + z * s v[j][i] = -x * s + z * c } } z = Math.sqrt(f * f + h * h) q[i - 1] = z c = f / z s = h / z f = c * g + s * y x = -s * g + c * y if (withu) { for (j = 0; j < m; j++) { y = u[j][i - 1] z = u[j][i] u[j][i - 1] = y * c + z * s u[j][i] = -y * s + z * c } } } e[l] = 0 e[k] = f q[k] = x } } // Number below eps should be zero for (i = 0; i < n; i++) { if (q[i] < eps) q[i] = 0 } return { u, q, v } }30:T1ece,/** SVD procedure as explained in "Singular Value Decomposition and Least Squares Solutions. By G.H. Golub et al." * * This procedure computes the singular values and complete orthogonal decomposition of a real rectangular matrix A: * A = U * diag(q) * V(t), U(t) * U = V(t) * V = I * where the arrays a, u, v, q represent A, U, V, q respectively. The actual parameters corresponding to a, u, v may * all be identical unless withu = withv = {true}. In this case, the actual parameters corresponding to u and v must * differ. m >= n is assumed (with m = a.length and n = a[0].length) * * @param a {Array} Represents the matrix A to be decomposed * @param [withu] {bool} {true} if U is desired {false} otherwise * @param [withv] {bool} {true} if U is desired {false} otherwise * @param [eps] {Number} A constant used in the test for convergence; should not be smaller than the machine precision * @param [tol] {Number} A machine dependent constant which should be set equal to B/eps0 where B is the smallest * positive number representable in the computer * * @returns {Object} An object containing: * q: A vector holding the singular values of A; they are non-negative but not necessarily ordered in * decreasing sequence * u: Represents the matrix U with orthonormalized columns (if withu is {true} otherwise u is used as * a working storage) * v: Represents the orthogonal matrix V (if withv is {true}, otherwise v is not used) * */ SVD = (a, withu, withv, eps, tol) => { // Define default parameters withu = withu !== undefined ? withu : true withv = withv !== undefined ? withv : true eps = eps || Math.pow(2, -52) tol = 1e-64 / eps // throw error if a is not defined if (!a) { throw new TypeError('Matrix a is not defined') } // Householder's reduction to bidiagonal form const n = a[0].length const m = a.length if (m < n) { throw new TypeError('Invalid matrix: m < n') } let i, j, k, l, l1, c, f, g, h, s, x, y, z g = 0 x = 0 const e = [] const u = [] const v = [] const mOrN = (withu === 'f') ? m : n // Initialize u for (i = 0; i < m; i++) { u[i] = new Array(mOrN).fill(0) } // Initialize v for (i = 0; i < n; i++) { v[i] = new Array(n).fill(0) } // Initialize q const q = new Array(n).fill(0) // Copy array a in u for (i = 0; i < m; i++) { for (j = 0; j < n; j++) { u[i][j] = a[i][j] } } for (i = 0; i < n; i++) { e[i] = g s = 0 l = i + 1 for (j = i; j < m; j++) { s += Math.pow(u[j][i], 2) } if (s < tol) { g = 0 } else { f = u[i][i] g = f < 0 ? Math.sqrt(s) : -Math.sqrt(s) h = f * g - s u[i][i] = f - g for (j = l; j < n; j++) { s = 0 for (k = i; k < m; k++) { s += u[k][i] * u[k][j] } f = s / h for (k = i; k < m; k++) { u[k][j] = u[k][j] + f * u[k][i] } } } q[i] = g s = 0 for (j = l; j < n; j++) { s += Math.pow(u[i][j], 2) } if (s < tol) { g = 0 } else { f = u[i][i + 1] g = f < 0 ? Math.sqrt(s) : -Math.sqrt(s) h = f * g - s u[i][i + 1] = f - g for (j = l; j < n; j++) { e[j] = u[i][j] / h } for (j = l; j < m; j++) { s = 0 for (k = l; k < n; k++) { s += u[j][k] * u[i][k] } for (k = l; k < n; k++) { u[j][k] = u[j][k] + s * e[k] } } } y = Math.abs(q[i]) + Math.abs(e[i]) if (y > x) { x = y } } // Accumulation of right-hand transformations if (withv) { for (i = n - 1; i >= 0; i--) { if (g !== 0) { h = u[i][i + 1] * g for (j = l; j < n; j++) { v[j][i] = u[i][j] / h } for (j = l; j < n; j++) { s = 0 for (k = l; k < n; k++) { s += u[i][k] * v[k][j] } for (k = l; k < n; k++) { v[k][j] = v[k][j] + s * v[k][i] } } } for (j = l; j < n; j++) { v[i][j] = 0 v[j][i] = 0 } v[i][i] = 1 g = e[i] l = i } } // Accumulation of left-hand transformations if (withu) { if (withu === 'f') { for (i = n; i < m; i++) { for (j = n; j < m; j++) { u[i][j] = 0 } u[i][i] = 1 } } for (i = n - 1; i >= 0; i--) { l = i + 1 g = q[i] for (j = l; j < mOrN; j++) { u[i][j] = 0 } if (g !== 0) { h = u[i][i] * g for (j = l; j < mOrN; j++) { s = 0 for (k = l; k < m; k++) { s += u[k][i] * u[k][j] } f = s / h for (k = i; k < m; k++) { u[k][j] = u[k][j] + f * u[k][i] } } for (j = i; j < m; j++) { u[j][i] = u[j][i] / g } } else { for (j = i; j < m; j++) { u[j][i] = 0 } } u[i][i] = u[i][i] + 1 } } // Diagonalization of the bidiagonal form eps = eps * x let testConvergence for (k = n - 1; k >= 0; k--) { for (let iteration = 0; iteration < 50; iteration++) { // test-f-splitting testConvergence = false for (l = k; l >= 0; l--) { if (Math.abs(e[l]) <= eps) { testConvergence = true break } if (Math.abs(q[l - 1]) <= eps) { break } } if (!testConvergence) { // cancellation of e[l] if l>0 c = 0 s = 1 l1 = l - 1 for (i = l; i < k + 1; i++) { f = s * e[i] e[i] = c * e[i] if (Math.abs(f) <= eps) { break // goto test-f-convergence } g = q[i] q[i] = Math.sqrt(f * f + g * g) h = q[i] c = g / h s = -f / h if (withu) { for (j = 0; j < m; j++) { y = u[j][l1] z = u[j][i] u[j][l1] = y * c + (z * s) u[j][i] = -y * s + (z * c) } } } } // test f convergence z = q[k] if (l === k) { // convergence if (z < 0) { // q[k] is made non-negative q[k] = -z if (withv) { for (j = 0; j < n; j++) { v[j][k] = -v[j][k] } } } break // break out of iteration loop and move on to next k value } // Shift from bottom 2x2 minor x = q[l] y = q[k - 1] g = e[k - 1] h = e[k] f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2 * h * y) g = Math.sqrt(f * f + 1) f = ((x - z) * (x + z) + h * (y / (f < 0 ? (f - g) : (f + g)) - h)) / x // Next QR transformation c = 1 s = 1 for (i = l + 1; i < k + 1; i++) { g = e[i] y = q[i] h = s * g g = c * g z = Math.sqrt(f * f + h * h) e[i - 1] = z c = f / z s = h / z f = x * c + g * s g = -x * s + g * c h = y * s y = y * c if (withv) { for (j = 0; j < n; j++) { x = v[j][i - 1] z = v[j][i] v[j][i - 1] = x * c + z * s v[j][i] = -x * s + z * c } } z = Math.sqrt(f * f + h * h) q[i - 1] = z c = f / z s = h / z f = c * g + s * y x = -s * g + c * y if (withu) { for (j = 0; j < m; j++) { y = u[j][i - 1] z = u[j][i] u[j][i - 1] = y * c + z * s u[j][i] = -y * s + z * c } } } e[l] = 0 e[k] = f q[k] = x } } // Number below eps should be zero for (i = 0; i < n; i++) { if (q[i] < eps) q[i] = 0 } return { u, q, v } }31:T60b,function svd(matrix, tol = 1e-10) { const A = matrix; const size = math.size(A); const m = size._data[0]; const n = size._data[1]; const U = math.zeros(m, m); const S = math.zeros(m, n); const V = math.zeros(n, n); // Step 1: Calculate A^T * A and A * A^T const ATA = math.multiply(math.transpose(A), A); const AAT = math.multiply(A, math.transpose(A)); // Step 2: Calculate the eigenvalues and eigenvectors of A^T * A const eigen = math.eigs(ATA); //Step 3: Sort eigenvalues and rearrange eigenvectors const sorted_indices = [...new Set(argSort(eigen.values.toArray()).reverse())]; const eigenvalues = eigen.values.subset(math.index(sorted_indices)); const eigenvectors = sorted_indices.map(i => eigen.eigenvectors[i].value); // Step 4: Calculate singular values and the pseudoinverse of S const singular_values = math.map(eigenvalues, math.sqrt); const non_zero_singular_values = math.map(singular_values, x => math.larger(x, tol)); const r = math.sum(non_zero_singular_values.map(x => x !== 0 ? 1 : 0)); // it seems to work up to here return math.range(1,r) S.subset( math.index(math.range(1, r), math.range(1, r)), math.diag(singular_values.subset(math.index(non_zero_singular_values))) ); return {eigenvalues, eigenvectors} const S_pseudo_inverse = math.pinv(S); // Step 5: Calculate U and V U.subset( [":", `:{r}`], math.multiply(A, eigenvectors.subset([":", `:{r}`])) ); V.subset([":", `:{r}`], eigenvectors.subset([":", `:{r}`])); return { U, S, V }; }32:T832,// Singular Value Decomposition (A = U Σ V^T) for math.js // Uses QR algorithm for eigenvalue decomposition of A^T A and AA^T // This is GROK code, doesn't work but might be a start. function svd2(A) { // Ensure input is a matrix const matrix = math.isMatrix(A) ? A : math.matrix(A); const size = matrix.size(); if (size.length !== 2) { throw new Error("Matrix must be two-dimensional"); } const m = size[0], n = size[1]; // Compute A^T A and AA^T const At = math.transpose(matrix); const AtA = math.multiply(At, matrix); // n × n const AAt = math.multiply(matrix, At); // m × m // Compute eigenvalues and eigenvectors using math.eigs // For V: eigenvectors of A^T A const eigsAtA = math.eigs(AtA); let singularValues = eigsAtA.values .toArray() .map((x) => Math.sqrt(Math.max(0, x))); // Ensure non-negative const V = eigsAtA.eigenvectors; // For U: eigenvectors of AA^T const eigsAAt = math.eigs(AAt); const U = eigsAAt.eigenvectors; // Sort singular values in descending order and adjust U, V const indices = singularValues .map((val, idx) => ({ val, idx })) .sort((a, b) => b.val - a.val) .map((item) => item.idx); singularValues = indices.map((i) => singularValues[i]); // Reorder V columns const V_array = V.isMatrix ? V.toArray() : V; const V_sorted = math.zeros(n, n); for (let j = 0; j < n; j++) { V_sorted.subset( math.index(math.range(0, n), j), math.subset(V, math.index(math.range(0, n), indices[j])) ); } // Reorder U columns const U_array = U.toArray(); const U_sorted = math.zeros(m, m); for (let j = 0; j < Math.min(m, n); j++) { U_sorted.subset( math.index(math.range(0, m), j), math.subset(U, math.index(math.range(0, m), indices[j])) ); } // Construct Σ (diagonal matrix) const Sigma = math.zeros(m, n); for (let i = 0; i < Math.min(m, n); i++) { Sigma.subset(math.index(i, i), singularValues[i]); } // Return U, Σ, V^T as math.js matrices return { U: U_sorted, S: Sigma, Vt: math.transpose(V_sorted) }; }33:[] 2d:{"title":"SVD","cells":[{"id":0,"value":"# SVD\n\nSingular Value Decomposition\n\n# todo\n\n- compare to https://github.com/danilosalvati/svd-js","pinned":false,"hidden":false,"mode":"md"},{"id":12,"value":"test = svd(math.matrix(a))","pinned":true,"hidden":false,"mode":"ojs"},{"id":123,"value":"```\nfrom scipy.linalg import svd\n>>> svd([[1,2],[10,8]])\n(array([[-0.15738271, -0.98753769],\n [-0.98753769, 0.15738271]]),\n array([12.96701928, 0.92542471]),\n array([[-0.77371363, -0.63353549],\n [ 0.63353549, -0.77371363]]))\n```\n","pinned":false,"hidden":false,"mode":"md"},{"id":47,"value":"SVDJS.SVD(a)","pinned":true,"hidden":false,"mode":"ojs"},{"id":190,"value":"SVD3(a)","pinned":true,"hidden":false,"mode":"ojs"},{"id":183,"value":"svd2(a)","pinned":true,"hidden":false,"mode":"ojs"},{"id":43,"value":"function argSort(a) {\n const b = a.map((x, i) => [i, x]);\n b.sort((a, b) => (a[1] > b[1] ? 1 : a[1] < b[1] ? -1 : 0));\n return b.map((x) => x[0]);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":189,"value":"$2f","pinned":true,"hidden":false,"mode":"ojs"},{"id":193,"value":"$30","pinned":true,"hidden":false,"mode":"ojs"},{"id":8,"value":"$31","pinned":true,"hidden":false,"mode":"ojs"},{"id":179,"value":"$32","pinned":true,"hidden":false,"mode":"ojs"},{"id":148,"value":"math.version","pinned":true,"hidden":false,"mode":"ojs"},{"id":3,"value":"math = require(\"mathjs@14.6.0\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":166,"value":"SVDJS = require(\"svd-js\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":168,"value":"a = [[1,2],[10,8]]","pinned":true,"hidden":false,"mode":"ojs"},{"id":170,"value":"result = SVDJS.SVD(a)","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W33","stdlib":"1","titling":"auto"}