Add PCA to statistics

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

Author: kylecorry31

src/main/kotlin/com/kylecorry/sol/math/statistics/PrincipalComponentAnalysis.kt

@@ -0,0 +1,35 @@
+package com.kylecorry.sol.math.statistics
+
+import com.kylecorry.sol.math.algebra.Matrix
+import com.kylecorry.sol.math.algebra.LinearAlgebra
+import com.kylecorry.sol.math.algebra.dot
+import com.kylecorry.sol.math.algebra.subtract
+import com.kylecorry.sol.math.algebra.transpose
+
+internal object PrincipalComponentAnalysis {
+
+    fun getComponents(x: Matrix, k: Int): Matrix {
+        require(k >= 0) { "k must be non-negative" }
+        require(k <= x.columns()) { "k must be less than or equal to the number of features" }
+
+        if (k == 0 || x.rows() == 0) {
+            return Matrix.zeros(x.rows(), k)
+        }
+
+        val featureMean = Matrix.create(1, x.columns()) { _, column ->
+            x.getColumn(column).average().toFloat()
+        }
+        val xTilde = x.subtract(featureMean)
+        val x2 = xTilde.transpose().dot(xTilde)
+        val eigen = LinearAlgebra.eigen(x2, TOLERANCE, MAX_ITERATIONS)
+        val sortedIndices = (0..<eigen.values.size).sortedByDescending { eigen.values[it] }
+        val components = Matrix.create(eigen.vectors.rows(), k) { row, column ->
+            eigen.vectors[row, sortedIndices[column]]
+        }
+        return xTilde.dot(components)
+    }
+
+    private const val TOLERANCE = 1e-6f
+    private const val MAX_ITERATIONS = 100
+
+}

src/main/kotlin/com/kylecorry/sol/math/statistics/Statistics.kt

@@ -433,4 +433,14 @@ object Statistics {
 
         return quantile(values, 0.75f, interpolate) - quantile(values, 0.25f, interpolate)
     }
+
+    /**
+     * Computes the PCA of the input matrix.
+     * @param x the input matrix where rows are samples and columns are features
+     * @param k the number of components to extract
+     * @return a matrix where rows are samples and columns are principal component projections (centered)
+     */
+    fun principalComponentAnalysis(x: Matrix, k: Int): Matrix {
+        return PrincipalComponentAnalysis.getComponents(x, k)
+    }
 }

src/test/kotlin/com/kylecorry/sol/math/statistics/StatisticsTest.kt

@@ -137,6 +137,50 @@ internal class StatisticsTest {
         assertEquals(expected, actual, 0.00001f)
     }
 
+    @Test
+    fun principalComponentAnalysis() {
+        val data = Matrix.create(
+            arrayOf(
+                floatArrayOf(2f, 1f, 5f),
+                floatArrayOf(4f, 3f, 4f),
+                floatArrayOf(6f, 2f, 3f),
+                floatArrayOf(8f, 5f, 2f),
+                floatArrayOf(10f, 7f, 1f)
+            )
+        )
+
+        val actual = Statistics.principalComponentAnalysis(data, 2)
+
+        assertEquals(5, actual.rows())
+        assertEquals(2, actual.columns())
+        assertProjectedEquals(
+            Matrix.create(
+                arrayOf(
+                    floatArrayOf(-5.164735f, 0.292427f),
+                    floatArrayOf(-2.196894f, 0.730517f),
+                    floatArrayOf(-0.881081f, -1.335551f),
+                    floatArrayOf(2.637435f, -0.062742f),
+                    floatArrayOf(5.605276f, 0.375348f)
+                )
+            ),
+            actual
+        )
+    }
+
+    private fun assertProjectedEquals(
+        expected: Matrix,
+        actual: Matrix,
+        tolerance: Float = 0.0001f
+    ) {
+        assertEquals(expected.rows(), actual.rows())
+        assertEquals(expected.columns(), actual.columns())
+        for (row in 0..<expected.rows()) {
+            for (col in 0..<expected.columns()) {
+                assertEquals(expected[row, col], actual[row, col], tolerance)
+            }
+        }
+    }
+
     companion object {
 
         @JvmStatic