Skip to content

Instantly share code, notes, and snippets.

@bellbind
Last active December 27, 2021 17:37
Show Gist options
  • Save bellbind/60fb876842346864baaf340cb789ad10 to your computer and use it in GitHub Desktop.
Save bellbind/60fb876842346864baaf340cb789ad10 to your computer and use it in GitHub Desktop.
[javascript] compute eigenvalues with QR method and Householder Trans
"use strict";
// helpers
function mat2d(mat) {
const n = Math.sqrt(mat.length);
return Array.from(Array(n), (_, i) => mat.slice(i * n, i * n + n));
}
function range(n, f = v => v) {
return Array.from(Array(n), (_, i) => f(i));
}
function sum(a, f = v => v) {
return a.reduce((r, v, i) => r + f(v, i), 0);
}
// returns hessenberg matrix: m[y][x] == 0 when y >= x + 2
function householder(mat) {
const n = Math.sqrt(mat.length);
const m = Array.from(mat);
console.assert(Number.isInteger(n));
const idx = (x, y) => y * n + x;
for (let k = 0; k < n - 2; k++) {
if (Math.abs(m[idx(k, k + 1)]) < Number.EPSILON) continue;
const u = range(n, i => i <= k ? 0 : m[idx(k, i)]);
const sigma = Math.sign(u[k + 1]) * Math.hypot(...u);
u[k + 1] += sigma;
const norm = Math.sqrt(2 * sigma * u[k + 1]);
const h = u.map(v => v / norm);
const sdx = range(n, i => sum(h, (v, j) => v * m[idx(j, i)]));
const sdy = range(n, i => sum(h, (v, j) => v * m[idx(i, j)]));
const hdx = sum(sdx, (v, i) => v * h[i]);
const hdy = sum(sdy, (v, i) => v * h[i]);
const dx = sdx.map((v, i) => 2 * (v - hdx * h[i]));
const dy = sdy.map((v, i) => 2 * (v - hdy * h[i]));
for (let y = 0; y < n; y++) for (let x = 0; x < n; x++) {
m[idx(x, y)] -= h[y] * dy[x] + h[x] * dx[y];
}
}
return m;
}
// returns upper triangler matrix: m[y][x] == 0 when y >= x + 1
function qr(mat) {
const n = Math.sqrt(mat.length);
const m = Array.from(mat);
console.assert(Number.isInteger(n));
const idx = (x, y) => y * n + x;
let k = n;
while (k >= 2) {
//console.log(k, m[idx(k - 2, k - 1)]);
//console.log(mat2d(m));
if (Math.abs(m[idx(k - 2, k - 1)]) < Number.EPSILON) {
k--;
continue;
}
// eigenvalue of last 2x2 sub matrix: l^2 - tr * l + det = 0
const a = m[idx(k - 1, k - 1)], b = m[idx(k - 2, k - 1)],
c = m[idx(k - 1, k - 2)], d = m[idx(k - 2, k - 2)];
const tr = a * d, det = a * d - b * c;
const disc = Math.sqrt(tr * tr - 4 * det) || 0;
const l1 = (tr + disc) / 2, l2 = (tr - disc) / 2;
const mu = a - (Math.abs(l1) < Math.abs(l2) ? l1 : l2);
// (option) pre process M = M - mu * I
for (let i = 0; i < k; i++) m[idx(i, i)] -= mu;
// init q as unit matrix
const idxq = (x, y) => y * k + x;
const q = Array(k * k).fill(0);
for (let i = 0; i < k; i++) q[idxq(i, i)] = 1;
// rotate to makes M => Q * R (R store to M)
for (let i = 0; i < k - 1; i++) {
const a1 = m[idx(i, i)], a2 = m[idx(i, i + 1)];
const base = Math.hypot(a1, a2);
const cos = base < Number.EPSILON ? 0 : a1 / base;
const sin = base < Number.EPSILON ? 0 : a2 / base;
// make R
m[idx(i, i)] = base;
m[idx(i, i + 1)] = 0;
for (let x = i + 1; x < k; x++) {
const e1 = m[idx(x, i)], e2 = m[idx(x, i + 1)];
m[idx(x, i)] = e1 * cos + e2 * sin;
m[idx(x, i + 1)] = e2 * cos - e1 * sin;
}
// make Q
for (let y = 0; y < k; y++) {
const e1 = q[idxq(i, y)], e2 = q[idxq(i + 1, y)];
q[idxq(i, y)] = e1 * cos + e2 * sin;
q[idxq(i + 1, y)] = e2 * cos - e1 * sin;
}
}
// next M as R * Q
for (let y = 0; y < k; y++) {
const ry = Array.from(Array(k - y), (_, j) => m[idx(y + j, y)]);
for (let x = 0; x < k; x++) {
m[idx(x, y)] = sum(ry, (v, j) => v * q[idxq(x, j + y)]);
}
}
// (option) post process M = M + mu * I
for (let i = 0; i < k; i++) m[idx(i, i)] += mu;
}
return m;
}
// list of eigen values square matrix (allow non symmetric)
function eigenvalues(mat) {
const ut = qr(householder(mat));
const n = Math.sqrt(ut.length);
return range(n, i => ut[i * n + i]);
}
// example:
// m = numpy.mat([[4, -6, 5], [-6, 3, 4], [5, 4, -3]])
// numpy.linalg.eigvals(m) #=> array([-9.12030391, 9.62192181, 3.4983821 ])
console.log(eigenvalues([
4, -6, 5,
-6, 3, 4,
5, 4, -3
]));
@bellbind
Copy link
Author

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment