/*
* Copyright (c) Huawei Technologies Co., Ltd. 2024-2024. All rights reserved.
*/
package matrix4cj
import std.math.*
/** Eigenvalues and eigenvectors of a real matrix.
<P>
If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is
diagonal and the eigenvector matrix V is orthogonal.
I.e. A = V.times(D.times(V.transpose())) and
V.times(V.transpose()) equals the identity matrix.
<P>
If A is not symmetric, then the eigenvalue matrix D is block diagonal
with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The
columns of V represent the eigenvectors in the sense that A*V = V*D,
i.e. A.times(V) equals V.times(D). The matrix V may be badly
conditioned, or even singular, so the validity of the equation
A = V*D*inverse(V) depends upon V.cond().
**/
public class EigenvalueDecomposition {
/**Row and column dimension (square matrix)*/
private var n: Int64 = 0
/**Symmetry flag*/
private var issymmetric: Bool = false
/**Arrays for internal storage of eigenvalues*/
private var d: Array<Float64>
/**Array for internal storage of eigenvectors*/
private var e: Array<Float64>
/**Array for internal storage of nonsymmetric Hessenberg form*/
private var V: Matrix
/**Working storage for nonsymmetric algorithm*/
private var H: Array<Array<Float64>> = Array<Array<Float64>>()
private var ort: Array<Float64> = Array<Float64>(0, repeat: 0.0)
private func tred2(): Unit {
// This is derived from the Algol procedures tred2 by
// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
// Fortran subroutine in EISPACK.
for (j in 0..n) {
d[j] = V[n - 1, j]
}
for (i in (n - 1)..0 : -1) {
// Scale to avoid under/overflow.
var scale: Float64 = 0.0
var h: Float64 = 0.0
for (k in 0..i) {
scale = scale + abs(d[k])
}
if (scale == 0.0) {
e[i] = d[i - 1]
for (j in 0..i) {
d[j] = V[i - 1, j]
V[i, j] = 0.0
V[j, i] = 0.0
}
} else {
// Generate Householder vector.
for (k in 0..i) {
d[k] /= scale
h += d[k] * d[k]
}
var f: Float64 = d[i - 1]
var g: Float64 = sqrt(h)
if (f > 0.0) {
g = -g
}
e[i] = scale * g
h = h - f * g
d[i - 1] = f - g
for (j in 0..i) {
e[j] = 0.0
}
// Apply similarity transformation to remaining columns.
for (j in 0..i) {
f = d[j]
V[j, i] = f
g = e[j] + V[j, j] * f
do {
var k = j + 1
while (k <= i - 1) {
g += V[k, j] * d[k]
e[k] += V[k, j] * f
k++
}
} while (false)
e[j] = g
}
f = 0.0
for (j in 0..i) {
e[j] /= h
f += e[j] * d[j]
}
let hh: Float64 = f / (h + h)
for (j in 0..i) {
e[j] -= hh * d[j]
}
for (j in 0..i) {
f = d[j]
g = e[j]
for (k in j..(i)) {
V[k, j] -= (f * e[k] + g * d[k])
}
d[j] = V[i - 1, j]
V[i, j] = 0.0
}
}
d[i] = h
}
// Accumulate transformations.
for (i in 0..(n - 1)) {
V[n - 1, i] = V[i, i]
V[i, i] = 1.0
let h = d[i + 1]
if (h != 0.0) {
for (k in 0..(i + 1)) {
d[k] = V[k, i + 1] / h
}
for (j in 0..i + 1) {
var g = 0.0
for (k in 0..(i + 1)) {
g += V[k, i + 1] * V[k, j]
}
for (k in 0..(i + 1)) {
V[k, j] -= g * d[k]
}
}
}
for (k in 0..(i + 1)) {
V[k, i + 1] = 0.0
}
}
for (j in 0..n) {
d[j] = V[n - 1, j]
V[n - 1, j] = 0.0
}
V[n - 1, n - 1] = 1.0
e[0] = 0.0
}
// Symmetric tridiagonal QL algorithm.
private func tql2(): Unit {
// This is derived from the Algol procedures tql2, by
// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
// Fortran subroutine in EISPACK.
for (i in 1..n) {
e[i - 1] = e[i]
}
e[n - 1] = 0.0
var f: Float64 = 0.0
var tst1: Float64 = 0.0
let eps: Float64 = pow(2.0, -52.0)
for (l in 0..n) {
tst1 = max(tst1, abs(d[Int64(l)]) + abs(e[Int64(l)]))
var m = l
while (m < n) {
if (abs(e[m]) <= eps * tst1) {
break
}
m++
}
if (m > l) {
var iter: Int64 = 0
do {
iter = iter + 1
var g: Float64 = d[l]
var p: Float64 = (d[l + 1] - g) / (2.0 * e[l])
var r: Float64 = Maths.hypot(p, 1.0)
if (p < 0.0) {
r = -r
}
d[l] = e[l] / (p + r)
d[l + 1] = e[l] * (p + r)
let dl1: Float64 = d[l + 1]
var h: Float64 = g - d[l]
for (i in (l + 2)..n) {
d[i] -= h
}
f = f + h
p = d[m]
var c: Float64 = 1.0
var c2: Float64 = c
var c3: Float64 = c
let el1: Float64 = e[l + 1]
var s: Float64 = 0.0
var s2: Float64 = 0.0
for (i in (m - 1)..l - 1 : -1) {
c3 = c2
c2 = c
s2 = s
g = c * e[i]
h = c * p
r = Maths.hypot(p, e[i])
e[i + 1] = s * r
s = e[i] / r
c = p / r
p = c * d[i] - s * g
d[i + 1] = h + s * (c * g + s * d[i])
for (k in 0..n) {
h = V[k, i + 1]
V[k, i + 1] = s * V[k, i] + c * h
V[k, i] = c * V[k, i] - s * h
}
}
p = -s * s2 * c3 * el1 * e[l] / dl1
e[l] = s * p
d[l] = c * p
} while (abs(e[l]) > eps * tst1)
}
d[l] = d[l] + f
e[l] = 0.0
}
// Sort eigenvalues and corresponding vectors.
for (i in 0..(n - 1)) {
var k = i
var p = d[i]
for (j in (i + 1)..n) {
if (d[j] < p) {
k = j
p = d[j]
}
}
if (k != i) {
d[k] = d[i]
d[i] = p
for (j in 0..n) {
p = V[j, i]
V[j, i] = V[j, k]
V[j, k] = p
}
}
}
}
// Nonsymmetric reduction to Hessenberg form.
private func orthes(): Unit {
// This is derived from the Algol procedures orthes and ortran,
// by Martin and Wilkinson, Handbook for Auto. Comp.,
// Vol.ii-Linear Algebra, and the corresponding
// Fortran subroutines in EISPACK.
let low: Int64 = 0
let high: Int64 = n - 1
for (m in low + 1..high) {
var scale = 0.0
for (i in m..high + 1) {
scale = scale + abs(H[i][m - 1])
}
if (scale != 0.0) {
var h: Float64 = 0.0
for (i in high..=m : -1) {
ort[i] = H[i][m - 1] / scale
h += ort[i] * ort[i]
}
var g: Float64 = sqrt(h)
if (ort[m] > 0.0) {
g = -g
}
h = h - ort[m] * g
ort[m] = ort[m] - g
for (j in m..n) {
var f: Float64 = 0.0
for (i in high..=m : -1) {
f += ort[i] * H[i][j]
}
f = f / h
for (i in m..(high + 1)) {
H[i][j] -= f * ort[i]
}
}
for (i in 0..(high + 1)) {
var f: Float64 = 0.0
for (j in high..=m : -1) {
f += ort[j] * H[i][j]
}
f = f / h
for (j in m..(high + 1)) {
H[i][j] -= f * ort[j]
}
}
ort[m] = scale * ort[m]
H[m][m - 1] = scale * g
}
}
for (i in 0..n) {
for (j in 0..n) {
V[i, j] = (if (i == j) {
1.0
} else {
0.0
})
}
}
for (m in (high - 1)..=(low + 1) : -1) {
if (H[m][m - 1] != 0.0) {
for (i in (m + 1)..(high + 1)) {
ort[i] = H[i][m - 1]
}
for (j in m..(high + 1)) {
var g: Float64 = 0.0
for (i in m..(high + 1)) {
g += ort[i] * V[i, j]
}
g = (g / ort[m]) / H[m][m - 1]
for (i in m..(high + 1)) {
V[i, j] += g * ort[i]
}
}
}
}
}
private var cdivr: Float64 = 0.0
private var cdivi: Float64 = 0.0
// Complex scalar division.
private func cdiv(xr: Float64, xi: Float64, yr: Float64, yi: Float64): Unit {
let r: Float64
let d: Float64
if (abs(yr) > abs(yi)) {
r = yi / yr
d = yr + r * yi
cdivr = (xr + r * xi) / d
cdivi = (xi - r * xr) / d
} else {
r = yr / yi
d = yi + r * yr
cdivr = (r * xr + xi) / d
cdivi = (r * xi - xr) / d
}
}
// Nonsymmetric reduction from Hessenberg to real Schur form.
private func hqr2(): Unit {
// This is derived from the Algol procedure hqr2,
// by Martin and Wilkinson, Handbook for Auto. Comp.,
// Vol.ii-Linear Algebra, and the corresponding
// Fortran subroutine in EISPACK.
let nn: Int64 = this.n
var n: Int64 = nn - 1
let low: Int64 = 0
let high: Int64 = nn - 1
let eps: Float64 = pow(2.0, -52.0)
var exshift: Float64 = 0.0
var p: Float64 = 0.0
var q: Float64 = 0.0
var r: Float64 = 0.0
var s: Float64 = 0.0
var z: Float64 = 0.0
var t: Float64 = 0.0
var w: Float64 = 0.0
var x: Float64 = 0.0
var y: Float64 = 0.0
// Store roots isolated by balanc and compute matrix norm
var norm: Float64 = 0.0
for (i in 0..nn) {
if (i < low || i > high) {
d[i] = H[i][i]
e[i] = 0.0
}
for (j in max(i - 1, 0)..nn) {
norm = norm + abs(H[i][j])
}
}
// Outer loop over eigenvalue index
var iter = 0
while (n >= low) {
// Look for single small sub-diagonal element
var l = n
while (l > low) {
s = abs(H[l - 1][l - 1]) + abs(H[l][l])
if (s == 0.0) {
s = norm
}
if (abs(H[l][l - 1]) < eps * s) {
break
}
l--
}
// Check for convergence
// One root found
if (l == n) {
H[n][n] = H[n][n] + exshift
d[n] = H[n][n]
e[n] = 0.0
n--
iter = 0
// Two roots found
} else if (l == n - 1) {
w = H[n][n - 1] * H[n - 1][n]
p = (H[n - 1][n - 1] - H[n][n]) / 2.0
q = p * p + w
z = sqrt(abs(q))
H[n][n] = H[n][n] + exshift
H[n - 1][n - 1] = H[n - 1][n - 1] + exshift
x = H[n][n]
// Real pair
if (q >= 0.0) {
if (p >= 0.0) {
z = p + z
} else {
z = p - z
}
d[n - 1] = x + z
d[n] = d[n - 1]
if (z != 0.0) {
d[n] = x - w / z
}
e[n - 1] = 0.0
e[n] = 0.0
x = H[n][n - 1]
s = abs(x) + abs(z)
p = x / s
q = z / s
r = sqrt(p * p + q * q)
p = p / r
q = q / r
// Row modification
for (j in (n - 1)..nn) {
z = H[n - 1][j]
H[n - 1][j] = q * z + p * H[n][j]
H[n][j] = q * H[n][j] - p * z
}
// Column modification
for (i in 0..(n + 1)) {
z = H[i][n - 1]
H[i][n - 1] = q * z + p * H[i][n]
H[i][n] = q * H[i][n] - p * z
}
// Accumulate transformations
for (i in low..(high + 1)) {
z = V[i, n - 1]
V[i, n - 1] = q * z + p * V[i, n]
V[i, n] = q * V[i, n] - p * z
}
// Complex pair
} else {
d[n - 1] = x + p
d[n] = x + p
e[n - 1] = z
e[n] = -z
}
n = n - 2
iter = 0
// No convergence yet
} else {
x = H[n][n]
y = 0.0
w = 0.0
if (l < n) {
y = H[n - 1][n - 1]
w = H[n][n - 1] * H[n - 1][n]
}
// Wilkinson's original ad hoc shift
if (iter == 10) {
exshift += x
for (i in low..(n + 1)) {
H[i][i] -= x
}
s = abs(H[n][n - 1]) + abs(H[n - 1][n - 2])
y = 0.75 * s
x = y
w = -0.4375 * s * s
}
// MATLAB's new ad hoc shift
if (iter == 30) {
s = (y - x) / 2.0
s = s * s + w
if (s > 0.0) {
s = sqrt(s)
if (y < x) {
s = -s
}
s = x - w / ((y - x) / 2.0 + s)
for (i in low..(n + 1)) {
H[i][i] -= s
}
exshift += s
w = 0.964
y = w
x = y
}
}
iter = iter + 1 // (Could check iteration count here.)
var m = n - 2
// Look for two consecutive small sub-diagonal elements
while (m >= l) {
z = H[m][m]
r = x - z
s = y - z
p = (r * s - w) / H[m + 1][m] + H[m][m + 1]
q = H[m + 1][m + 1] - z - r - s
r = H[m + 2][m + 1]
s = abs(p) + abs(q) + abs(r)
p = p / s
q = q / s
r = r / s
if (m == l) {
break
}
if (abs(H[m][m - 1]) * (abs(q) + abs(r)) < eps * (abs(p) * (abs(H[m - 1][m - 1]) + abs(z) + abs(
H[m + 1][m + 1])))) {
break
}
m--
}
for (i in (m + 2)..(n + 1)) {
H[i][i - 2] = 0.0
if (i > m + 2) {
H[i][i - 3] = 0.0
}
}
// Double QR step involving rows l:n and columns m:n
for (k in m..n) {
let notlast = (k != n - 1)
if (k != m) {
p = H[k][k - 1]
q = H[k + 1][k - 1]
r = (if (notlast) {
H[k + 2][k - 1]
} else {
0.0
})
x = abs(p) + abs(q) + abs(r)
if (x == 0.0) {
continue
}
p = p / x
q = q / x
r = r / x
}
s = sqrt(p * p + q * q + r * r)
if (p < 0.0) {
s = -s
}
if (s != 0.0) {
if (k != m) {
H[k][k - 1] = -s * x
} else if (l != m) {
H[k][k - 1] = -H[k][k - 1]
}
p = p + s
x = p / s
y = q / s
z = r / s
q = q / p
r = r / p
// Row modification
for (j in k..nn) {
p = H[k][j] + q * H[k + 1][j]
if (notlast) {
p = p + r * H[k + 2][j]
H[k + 2][j] = H[k + 2][j] - p * z
}
H[k][j] = H[k][j] - p * x
H[k + 1][j] = H[k + 1][j] - p * y
}
// Column modification
for (i in 0..=min(n, k + 3)) {
p = x * H[i][k] + y * H[i][k + 1]
if (notlast) {
p = p + z * H[i][k + 2]
H[i][k + 2] = H[i][k + 2] - p * r
}
H[i][k] = H[i][k] - p
H[i][k + 1] = H[i][k + 1] - p * q
}
// Accumulate transformations
for (i in low..(high + 1)) {
p = x * V[i, k] + y * V[i, k + 1]
if (notlast) {
p = p + z * V[i, k + 2]
V[i, k + 2] = V[i, k + 2] - p * r
}
V[i, k] = V[i, k] - p
V[i, k + 1] = V[i, k + 1] - p * q
}
}
}
}
}
// Backsubstitute to find vectors of upper triangular form
if (norm == 0.0) {
return
}
for (n_t in nn - 1..=0 : -1) {
n = n_t
p = d[n]
q = e[n]
if (q == 0.0) {
var l: Int64 = n
H[n][n] = 1.0
for (i in (n - 1)..-1 : -1) {
w = H[i][i] - p
r = 0.0
for (j in l..(n + 1)) {
r = r + H[i][j] * H[j][n]
}
if (e[i] < 0.0) {
z = w
s = r
} else {
l = i
if (e[Int64(i)] == 0.0) {
if (w != 0.0) {
H[i][n] = -r / w
} else {
H[i][n] = -r / (eps * norm)
}
// Solve real equations
} else {
x = H[i][i + 1]
y = H[i + 1][i]
q = (d[i] - p) * (d[i] - p) + e[i] * e[i]
t = (x * s - z * r) / q
H[i][n] = t
if (abs(x) > abs(z)) {
H[i + 1][n] = (-r - w * t) / x
} else {
H[i + 1][n] = (-s - y * t) / z
}
}
// Overflow control
t = abs(H[i][n])
if ((eps * t) * t > 1.0) {
for (j in i..(n + 1)) {
H[j][n] = H[j][n] / t
}
}
}
}
// Complex vector
} else if (q < 0.0) {
var l: Int64 = n - 1
// Last vector component imaginary so matrix is triangular
if (abs(H[n][n - 1]) > abs(H[n - 1][n])) {
H[n - 1][n - 1] = q / H[n][n - 1]
H[n - 1][n] = -(H[n][n] - p) / H[n][n - 1]
} else {
cdiv(0.0, -H[n - 1][n], H[n - 1][n - 1] - p, q)
H[n - 1][n - 1] = cdivr
H[n - 1][n] = cdivi
}
H[n][n - 1] = 0.0
H[n][n] = 1.0
for (i in (n - 2)..=0 : -1) {
var ra: Float64 = 0.0
var sa: Float64 = 0.0
var vr: Float64 = 0.0
let vi: Float64
ra = 0.0
sa = 0.0
do {
var j = l
while (j <= n) {
ra = ra + H[i][j] * H[j][n - 1]
sa = sa + H[i][j] * H[j][n]
j++
}
} while (false)
w = H[i][i] - p
if (e[i] < 0.0) {
z = w
r = ra
s = sa
} else {
l = i
if (e[i] == 0.0) {
cdiv(-ra, -sa, w, q)
H[i][n - 1] = cdivr
H[i][n] = cdivi
} else {
// Solve complex equations
x = H[i][i + 1]
y = H[i + 1][i]
vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q
vi = (d[i] - p) * 2.0 * q
if (vr == 0.0 && vi == 0.0) {
vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(z))
}
cdiv(x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi)
H[i][n - 1] = cdivr
H[i][n] = cdivi
if (abs(x) > (abs(z) + abs(q))) {
H[i + 1][n - 1] = (-ra - w * H[i][n - 1] + q * H[i][n]) / x
H[i + 1][n] = (-sa - w * H[i][n] - q * H[i][n - 1]) / x
} else {
cdiv(-r - y * H[i][n - 1], -s - y * H[i][n], z, q)
H[i + 1][n - 1] = cdivr
H[i + 1][n] = cdivi
}
}
// Overflow control
t = max(abs(H[i][n - 1]), abs(H[i][n]))
if ((eps * t) * t > 1.0) {
for (j in i..(n + 1)) {
H[j][n - 1] = H[j][n - 1] / t
H[j][n] = H[j][n] / t
}
}
}
}
}
}
// Vectors of isolated roots
for (i in 0..nn) {
if (i < low || i > high) {
for (j in i..nn) {
V[i, j] = H[i][j]
}
}
}
// Back transformation to get eigenvectors of original matrix
for (j in (nn - 1)..(low - 1) : -1) {
for (i in low..(high + 1)) {
z = 0.0
for (k in low..=min(j, high)) {
z = z + V[i, k] * H[k][j]
}
V[i, j] = z
}
}
}
/**
* Check for symmetry, then construct the eigenvalue decomposition
* Structure to access D and V.
* @param Arg Square matrix
*/
public init(Arg: Matrix) {
n = Arg.colNum
V = Matrix(n, n)
d = Array<Float64>(n) {_ => 0.0}
e = Array<Float64>(n) {_ => 0.0}
issymmetric = true
for (j in 0..n where issymmetric) {
for (i in 0..n where issymmetric) {
issymmetric = (Arg[i, j] == Arg[j, i])
}
}
if (issymmetric) {
for (i in 0..n) {
for (j in 0..n) {
V[i, j] = Arg[i, j]
}
}
// Tridiagonalize.
tred2()
// Diagonalize.
tql2()
} else {
let array_dim_generated_n = n
H = Array<Array<Float64>>(array_dim_generated_n) {_ => Array<Float64>(array_dim_generated_n) {_ => 0.0}}
ort = Array<Float64>(n) {_ => 0.0}
for (j in 0..n) {
for (i in 0..n) {
H[i][j] = Arg[i, j]
}
}
// Reduce to Hessenberg form.
orthes()
// Reduce Hessenberg to real Schur form.
hqr2()
}
}
/**
* Return the eigenvector matrix
* @return V
*/
public func getV(): Matrix {
return V
}
/**
* Return the real parts of the eigenvalues
* @return real(diag(D))
*/
public func getRealEigenvalues(): Array<Float64> {
return d
}
/**
* Return the imaginary parts of the eigenvalues
* @return imag(diag(D))
*/
public func getImagEigenvalues(): Array<Float64> {
return e
}
/**
* Return the block diagonal eigenvalue matrix D
* @return D
*/
public func getD(): Matrix {
let X = Matrix(n, n, value: 0.0)
for (i in 0..n){
X[i, i] = d[i]
if (e[i] > 0.0) {
X[i, i + 1] = e[i]
} else if (e[i] < 0.0) {
X[i, i - 1] = e[i]
}
}
return X
}
}
