Calculate eigen values and vectors

Previously was only calculating values. Test made using numpy's eigen method.

Signed-off-by: Kyle Corry <kylecorry31@gmail.com>

Author: kylecorry31

src/main/kotlin/com/kylecorry/sol/math/algebra/LinearAlgebra.kt

@@ -264,7 +264,11 @@ object LinearAlgebra {
             }
 
             r[j, j] = norm(v)
-            q.setColumn(j, v.divide(r[j, j]).getRow(0))
+            if (Arithmetic.isZero(r[j, j])) {
+                q.setColumn(j, orthogonalUnitVector(q, j))
+            } else {
+                q.setColumn(j, v.divide(r[j, j]).getRow(0))
+            }
         }
 
         return q to r
@@ -279,21 +283,32 @@ object LinearAlgebra {
         }
     }
 
-    fun eigenvalues(m: Matrix, tolerance: Float = 1e-12f, maxIterations: Int = 1000): FloatArray {
+    fun eigen(
+        m: Matrix,
+        tolerance: Float = 1e-12f,
+        maxIterations: Int = 1000
+    ): EigenDecomposition {
+        require(m.rows() == m.columns()) { "Matrix must be square" }
+        require(m.rows() > 0) { "Matrix must not be empty" }
+        require(tolerance > 0f) { "Tolerance must be greater than zero" }
+        require(maxIterations >= 0) { "Max iterations must be non-negative" }
+
         var old = m
         var new = m
+        var vectors = Matrix.identity(m.rows())
 
         var diff = Float.MAX_VALUE
         var iterations = 0
         while (diff > tolerance && iterations < maxIterations) {
             val (q, r) = qr(new)
             new = r.dot(q)
+            vectors = vectors.dot(q)
             diff = max(new.subtract(old).mapped { it.absoluteValue })
             old = new
             iterations++
         }
 
-        return diagonal(new).getRow(0)
+        return EigenDecomposition(Vector(diagonal(new).getRow(0)), normalizeEigenvectorSigns(vectors))
     }
 
     fun norm(m: Matrix): Float {
@@ -451,6 +466,45 @@ object LinearAlgebra {
         return dot(a, x).toVector() - b
     }
 
+    private fun orthogonalUnitVector(q: Matrix, column: Int): FloatArray {
+        for (attempt in 0..<q.rows()) {
+            val vector = FloatArray(q.rows()) { index ->
+                if (index == (column + attempt) % q.rows()) 1f else 0f
+            }
+
+            for (previous in 0..<column) {
+                var projection = 0f
+                for (row in 0..<q.rows()) {
+                    projection += vector[row] * q[row, previous]
+                }
+
+                for (row in 0..<q.rows()) {
+                    vector[row] -= projection * q[row, previous]
+                }
+            }
+
+            val norm = sqrt(vector.sumOf { it * it.toDouble() }).toFloat()
+            if (!Arithmetic.isZero(norm)) {
+                return FloatArray(vector.size) { index -> vector[index] / norm }
+            }
+        }
+
+        return FloatArray(q.rows()) { index -> if (index == column % q.rows()) 1f else 0f }
+    }
+
+    private fun normalizeEigenvectorSigns(vectors: Matrix): Matrix {
+        return Matrix.create(vectors.rows(), vectors.columns()) { row, column ->
+            var maxIndex = 0
+            for (i in 1..<vectors.rows()) {
+                if (vectors[i, column].absoluteValue > vectors[maxIndex, column].absoluteValue) {
+                    maxIndex = i
+                }
+            }
+
+            vectors[row, column] * if (vectors[maxIndex, column] >= 0f) 1f else -1f
+        }
+    }
+
     fun appendColumn(m: Matrix, value: Float): Matrix {
         return Matrix.create(m.rows(), m.columns() + 1) { r, c ->
             if (c < m.columns()) {
@@ -482,4 +536,6 @@ object LinearAlgebra {
         }
     }
 
+    data class EigenDecomposition(val values: Vector, val vectors: Matrix)
+
 }

src/main/kotlin/com/kylecorry/sol/math/algebra/Matrix.kt

@@ -239,10 +239,18 @@ fun Matrix.add(other: Matrix): Matrix {
     return LinearAlgebra.add(this, other)
 }
 
+fun Matrix.add(value: Float): Matrix {
+    return LinearAlgebra.add(this, value)
+}
+
 fun Matrix.subtract(other: Matrix): Matrix {
     return LinearAlgebra.subtract(this, other)
 }
 
+fun Matrix.subtract(value: Float): Matrix {
+    return LinearAlgebra.subtract(this, value)
+}
+
 fun Matrix.multiply(other: Matrix): Matrix {
     return LinearAlgebra.multiply(this, other)
 }

src/test/kotlin/com/kylecorry/sol/math/algebra/LinearAlgebraTest.kt

@@ -505,6 +505,34 @@ class LinearAlgebraTest {
         assertEquals(sqrt(91.0).toFloat(), result, 0.00001f)
     }
 
+    @Test
+    fun eigen() {
+        val matrix = Matrix.create(
+            arrayOf(
+                floatArrayOf(4f, 1f, 1f),
+                floatArrayOf(1f, 3f, 0.5f),
+                floatArrayOf(1f, 0.5f, 2f)
+            )
+        )
+
+        val result = LinearAlgebra.eigen(matrix)
+
+        assertEquals(5.028030f, result.values[0], 0.0001f)
+        assertEquals(2.392581f, result.values[1], 0.0001f)
+        assertEquals(1.579389f, result.values[2], 0.0001f)
+        assertEquals(3, result.vectors.rows())
+        assertEquals(3, result.vectors.columns())
+        assertEquals(0.805083f, result.vectors[0, 0], 0.0001f)
+        assertEquals(0.482157f, result.vectors[1, 0], 0.0001f)
+        assertEquals(0.345493f, result.vectors[2, 0], 0.0001f)
+        assertEquals(-0.477410f, result.vectors[0, 1], 0.0001f)
+        assertEquals(0.872386f, result.vectors[1, 1], 0.0001f)
+        assertEquals(-0.104989f, result.vectors[2, 1], 0.0001f)
+        assertEquals(-0.352024f, result.vectors[0, 2], 0.0001f)
+        assertEquals(-0.080417f, result.vectors[1, 2], 0.0001f)
+        assertEquals(0.932530f, result.vectors[2, 2], 0.0001f)
+    }
+
     @Test
     fun max() {
         val m1 = Matrix.create(3, 2, 0f)