benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/SVDBenchmark.kt

@@ -8,7 +8,7 @@ import kotlinx.benchmark.Benchmark
 import kotlinx.benchmark.Blackhole
 import kotlinx.benchmark.Scope
 import kotlinx.benchmark.State
-import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.svd
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.svdGolubKahan
 import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
 import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.svdPowerMethod
 
@@ -26,9 +26,9 @@ class SVDBenchmark {
     }
 
     @Benchmark
-    fun svdGolabKahan(blackhole: Blackhole) {
+    fun svdGolubKahan(blackhole: Blackhole) {
         blackhole.consume(
-            tensor.svd()
+            tensor.svdGolubKahan()
         )
     }
 }

kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt

@@ -834,6 +834,10 @@ public open class DoubleTensorAlgebra :
     }
 
     override fun StructureND<Double>.svd(): Triple<DoubleTensor, DoubleTensor, DoubleTensor> {
+        return this.svdGolubKahan()
+    }
+
+    public fun StructureND<Double>.svdGolubKahan(): Triple<DoubleTensor, DoubleTensor, DoubleTensor> {
         val size = tensor.dimension
         val commonShape = tensor.shape.sliceArray(0 until size - 2)
         val (n, m) = tensor.shape.sliceArray(size - 2 until size)
@@ -855,7 +859,7 @@ public open class DoubleTensorAlgebra :
                     .slice(matrix.bufferStart until matrix.bufferStart + matrixSize)
                     .toDoubleArray()
             )
-            curMatrix.as2D().svdHelper(uTensors[index].as2D(), sTensorVectors[index], vTensors[index].as2D())
+            curMatrix.as2D().svdGolubKahanHelper(uTensors[index].as2D(), sTensorVectors[index], vTensors[index].as2D())
         }
 
         return Triple(uTensor.transpose(), sTensor, vTensor)

kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt

@@ -372,7 +372,7 @@ private fun SIGN(a: Double, b: Double): Double {
     else
         return -abs(a)
 }
-internal fun MutableStructure2D<Double>.svdHelper(u: MutableStructure2D<Double>, w: BufferedTensor<Double>, v: MutableStructure2D<Double>) {
+internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2D<Double>, w: BufferedTensor<Double>, v: MutableStructure2D<Double>) {
     val shape = this.shape
     val m = shape.component1()
     val n = shape.component2()

kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt

@@ -179,7 +179,7 @@ internal class TestDoubleLinearOpsTensorAlgebra {
             2.000000, 3.000000, 4.000000,
             3.000000, 4.000000, 5.000000,
             4.000000, 5.000000, 6.000000,
-            5.000000, 6.000000, 7.000000
+            5.000000, 6.000000, 9.000000
         )
         testSVDFor(fromArray(intArrayOf(5, 3), buffer1))
         val buffer2 = doubleArrayOf(
@@ -198,10 +198,52 @@ internal class TestDoubleLinearOpsTensorAlgebra {
         testSVDFor(tensor2)
     }
 
+    @Test
+    fun testSVDGolubKahan() = DoubleTensorAlgebra{
+        testSVDGolubKahanFor(fromArray(intArrayOf(2, 3), doubleArrayOf(1.0, 2.0, 3.0, 4.0, 5.0, 6.0)))
+        testSVDGolubKahanFor(fromArray(intArrayOf(2, 2), doubleArrayOf(-1.0, 0.0, 239.0, 238.0)))
+        val buffer1 = doubleArrayOf(
+            1.000000, 2.000000, 3.000000,
+            2.000000, 3.000000, 4.000000,
+            3.000000, 4.000000, 5.000000,
+            4.000000, 5.000000, 6.000000,
+            5.000000, 6.000000, 9.000000
+        )
+        testSVDGolubKahanFor(fromArray(intArrayOf(5, 3), buffer1))
+        val buffer2 = doubleArrayOf(
+            1.0, 2.0, 3.0, 2.0, 3.0,
+            4.0, 3.0, 4.0, 5.0, 4.0,
+            5.0, 6.0, 5.0, 6.0, 7.0
+        )
+        testSVDGolubKahanFor(fromArray(intArrayOf(3, 5), buffer2))
+    }
+
+    @Test
+    fun testBatchedSVDGolubKahan() = DoubleTensorAlgebra{
+        val tensor1 = randomNormal(intArrayOf(2, 5, 3), 0)
+        testSVDGolubKahanFor(tensor1)
+        val tensor2 = DoubleTensorAlgebra.randomNormal(intArrayOf(30, 30, 30), 0)
+        testSVDGolubKahanFor(tensor2)
+    }
+
     @Test
     fun testSVDPowerMethod() = DoubleTensorAlgebra{
         testSVDPowerMethodFor(fromArray(intArrayOf(2, 3), doubleArrayOf(1.0, 2.0, 3.0, 4.0, 5.0, 6.0)))
         testSVDPowerMethodFor(fromArray(intArrayOf(2, 2), doubleArrayOf(-1.0, 0.0, 239.0, 238.0)))
+        val buffer1 = doubleArrayOf(
+            1.000000, 2.000000, 3.000000,
+            2.000000, 3.000000, 4.000000,
+            3.000000, 4.000000, 5.000000,
+            4.000000, 5.000000, 6.000000,
+            5.000000, 6.000000, 9.000000
+        )
+        testSVDPowerMethodFor(fromArray(intArrayOf(5, 3), buffer1))
+        val buffer2 = doubleArrayOf(
+            1.0, 2.0, 3.0, 2.0, 3.0,
+            4.0, 3.0, 4.0, 5.0, 4.0,
+            5.0, 6.0, 5.0, 6.0, 7.0
+        )
+        testSVDPowerMethodFor(fromArray(intArrayOf(3, 5), buffer2))
     }
 
     @Test
@@ -237,6 +279,18 @@ private fun DoubleTensorAlgebra.testSVDFor(tensor: DoubleTensor) {
     assertTrue(tensor.eq(tensorSVD))
 }
 
+private fun DoubleTensorAlgebra.testSVDGolubKahanFor(tensor: DoubleTensor) {
+    val svd = tensor.svdGolubKahan()
+
+    val tensorSVD = svd.first
+        .dot(
+            diagonalEmbedding(svd.second)
+                .dot(svd.third.transpose())
+        )
+
+    assertTrue(tensor.eq(tensorSVD))
+}
+
 private fun DoubleTensorAlgebra.testSVDPowerMethodFor(tensor: DoubleTensor, epsilon: Double = 1e-10) {
     val svd = tensor.svdPowerMethod()