added new svd algorithm (Golub Kahan) and used by default for svd

kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt

@@ -19,9 +19,9 @@ import org.ejml.sparse.csc.factory.DecompositionFactory_DSCC
 import org.ejml.sparse.csc.factory.DecompositionFactory_FSCC
 import org.ejml.sparse.csc.factory.LinearSolverFactory_DSCC
 import org.ejml.sparse.csc.factory.LinearSolverFactory_FSCC
-import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.linear.*
 import space.kscience.kmath.linear.Matrix
+import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.nd.StructureFeature
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.FloatField

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

@@ -706,7 +706,7 @@ public open class DoubleTensorAlgebra :
     override fun svd(
         structureND: StructureND<Double>,
     ): Triple<StructureND<Double>, StructureND<Double>, StructureND<Double>> =
-        svd(structureND = structureND, epsilon = 1e-10)
+        svdGolubKahan(structureND = structureND, epsilon = 1e-10)
 
     override fun symEig(structureND: StructureND<Double>): Pair<DoubleTensor, DoubleTensor> =
         symEigJacobi(structureND = structureND, maxIteration = 50, epsilon = 1e-15)

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

@@ -7,10 +7,7 @@ package space.kscience.kmath.tensors.core.internal
 
 import space.kscience.kmath.nd.*
 import space.kscience.kmath.operations.invoke
-import space.kscience.kmath.structures.DoubleBuffer
+import space.kscience.kmath.structures.*
-import space.kscience.kmath.structures.IntBuffer
-import space.kscience.kmath.structures.asBuffer
-import space.kscience.kmath.structures.indices
 import space.kscience.kmath.tensors.core.*
 import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.div
 import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.dot
@@ -316,6 +313,302 @@ internal fun DoubleTensorAlgebra.svdHelper(
     }
 }
 
+private fun pythag(a: Double, b: Double): Double {
+    val at: Double = abs(a)
+    val bt: Double = abs(b)
+    val ct: Double
+    val result: Double
+    if (at > bt) {
+        ct = bt / at
+        result = at * sqrt(1.0 + ct * ct)
+    } else if (bt > 0.0) {
+        ct = at / bt
+        result = bt * sqrt(1.0 + ct * ct)
+    } else result = 0.0
+    return result
+}
+
+private fun SIGN(a: Double, b: Double): Double {
+    if (b >= 0.0)
+        return abs(a)
+    else
+        return -abs(a)
+}
+
+internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2D<Double>, w: BufferedTensor<Double>,
+                                                            v: MutableStructure2D<Double>, iterations: Int, epsilon: Double) {
+    val shape = this.shape
+    val m = shape.component1()
+    val n = shape.component2()
+    var f = 0.0
+    val rv1 = DoubleArray(n)
+    var s = 0.0
+    var scale = 0.0
+    var anorm = 0.0
+    var g = 0.0
+    var l = 0
+
+    val wStart = 0
+    val wBuffer = w.source
+
+    for (i in 0 until n) {
+        /* left-hand reduction */
+        l = i + 1
+        rv1[i] = scale * g
+        g = 0.0
+        s = 0.0
+        scale = 0.0
+        if (i < m) {
+            for (k in i until m) {
+                scale += abs(this[k, i]);
+            }
+            if (abs(scale) > epsilon) {
+                for (k in i until m) {
+                    this[k, i] = (this[k, i] / scale)
+                    s += this[k, i] * this[k, i]
+                }
+                f = this[i, i]
+                if (f >= 0) {
+                    g = (-1) * abs(sqrt(s))
+                } else {
+                    g = abs(sqrt(s))
+                }
+                val h = f * g - s
+                this[i, i] = f - g
+                if (i != n - 1) {
+                    for (j in l until n) {
+                        s = 0.0
+                        for (k in i until m) {
+                            s += this[k, i] * this[k, j]
+                        }
+                        f = s / h
+                        for (k in i until m) {
+                            this[k, j] += f * this[k, i]
+                        }
+                    }
+                }
+                for (k in i until m) {
+                    this[k, i] = this[k, i] * scale
+                }
+            }
+        }
+
+        wBuffer[wStart + i] = scale * g
+        /* right-hand reduction */
+        g = 0.0
+        s = 0.0
+        scale = 0.0
+        if (i < m && i != n - 1) {
+            for (k in l until n) {
+                scale += abs(this[i, k])
+            }
+            if (abs(scale) > epsilon) {
+                for (k in l until n) {
+                    this[i, k] = this[i, k] / scale
+                    s += this[i, k] * this[i, k]
+                }
+                f = this[i, l]
+                if (f >= 0) {
+                    g = (-1) * abs(sqrt(s))
+                } else {
+                    g = abs(sqrt(s))
+                }
+                val h = f * g - s
+                this[i, l] = f - g
+                for (k in l until n) {
+                    rv1[k] = this[i, k] / h
+                }
+                if (i != m - 1) {
+                    for (j in l until m) {
+                        s = 0.0
+                        for (k in l until n) {
+                            s += this[j, k] * this[i, k]
+                        }
+                        for (k in l until n) {
+                            this[j, k] += s * rv1[k]
+                        }
+                    }
+                }
+                for (k in l until n) {
+                    this[i, k] = this[i, k] * scale
+                }
+            }
+        }
+        anorm = max(anorm, (abs(wBuffer[wStart + i]) + abs(rv1[i])));
+    }
+
+    for (i in n - 1 downTo 0) {
+        if (i < n - 1) {
+            if (abs(g) > epsilon) {
+                for (j in l until n) {
+                    v[j, i] = (this[i, j] / this[i, l]) / g
+                }
+                for (j in l until n) {
+                    s = 0.0
+                    for (k in l until n)
+                        s += this[i, k] * v[k, j]
+                    for (k in l until n)
+                        v[k, j] += s * v[k, i]
+                }
+            }
+            for (j in l until n) {
+                v[i, j] = 0.0
+                v[j, i] = 0.0
+            }
+        }
+        v[i, i] = 1.0
+        g = rv1[i]
+        l = i
+    }
+
+    for (i in min(n, m) - 1 downTo 0) {
+        l = i + 1
+        g =  wBuffer[wStart + i]
+        for (j in l until n) {
+            this[i, j] = 0.0
+        }
+        if (abs(g) > epsilon) {
+            g = 1.0 / g
+            for (j in l until n) {
+                s = 0.0
+                for (k in l until m) {
+                    s += this[k, i] * this[k, j]
+                }
+                f = (s / this[i, i]) * g
+                for (k in i until m) {
+                    this[k, j] += f * this[k, i]
+                }
+            }
+            for (j in i until m) {
+                this[j, i] *= g
+            }
+        } else {
+            for (j in i until m) {
+                this[j, i] = 0.0
+            }
+        }
+        this[i, i] += 1.0
+    }
+
+    var flag = 0
+    var nm = 0
+    var c = 0.0
+    var h = 0.0
+    var y = 0.0
+    var z = 0.0
+    var x = 0.0
+    for (k in n - 1 downTo 0) {
+        for (its in 1 until iterations) {
+            flag = 1
+            for (newl in k downTo 0) {
+                nm = newl - 1
+                if (abs(rv1[newl]) + anorm == anorm) {
+                    flag = 0
+                    l = newl
+                    break
+                }
+                if (abs(wBuffer[wStart + nm]) + anorm == anorm) {
+                    l = newl
+                    break
+                }
+            }
+
+            if (flag != 0) {
+                c = 0.0
+                s = 1.0
+                for (i in l until k + 1) {
+                    f = s * rv1[i]
+                    rv1[i] = c * rv1[i]
+                    if (abs(f) + anorm == anorm) {
+                        break
+                    }
+                    g = wBuffer[wStart + i]
+                    h = pythag(f, g)
+                    wBuffer[wStart + i] = h
+                    h = 1.0 / h
+                    c = g * h
+                    s = (-f) * h
+                    for (j in 0 until m) {
+                        y = this[j, nm]
+                        z = this[j, i]
+                        this[j, nm] = y * c + z * s
+                        this[j, i] = z * c - y * s
+                    }
+                }
+            }
+
+            z =  wBuffer[wStart + k]
+            if (l == k) {
+                if (z < 0.0) {
+                    wBuffer[wStart + k] = -z
+                    for (j in 0 until n)
+                        v[j, k] = -v[j, k]
+                }
+                break
+            }
+
+            x =  wBuffer[wStart + l]
+            nm = k - 1
+            y =  wBuffer[wStart + nm]
+            g = rv1[nm]
+            h = rv1[k]
+            f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y)
+            g = pythag(f, 1.0)
+            f = ((x - z) * (x + z) + h * ((y / (f + SIGN(g, f))) - h)) / x
+            c = 1.0
+            s = 1.0
+
+            var i = 0
+            for (j in l until nm + 1) {
+                i = j + 1
+                g = rv1[i]
+                y = wBuffer[wStart + i]
+                h = s * g
+                g = c * g
+                z = pythag(f, h)
+                rv1[j] = z
+                c = f / z
+                s = h / z
+                f = x * c + g * s
+                g = g * c - x * s
+                h = y * s
+                y *= c
+
+                for (jj in 0 until n) {
+                    x = v[jj, j];
+                    z = v[jj, i];
+                    v[jj, j] = x * c + z * s;
+                    v[jj, i] = z * c - x * s;
+                }
+                z = pythag(f, h)
+                wBuffer[wStart + j] = z
+                if (abs(z) > epsilon) {
+                    z = 1.0 / z
+                    c = f * z
+                    s = h * z
+                }
+                f = c * g + s * y
+                x = c * y - s * g
+                for (jj in 0 until m) {
+                    y = this[jj, j]
+                    z = this[jj, i]
+                    this[jj, j] = y * c + z * s
+                    this[jj, i] = z * c - y * s
+                }
+            }
+            rv1[l] = 0.0
+            rv1[k] = f
+            wBuffer[wStart + k] = x
+        }
+    }
+
+    for (i in 0 until n) {
+        for (j in 0 until m) {
+            u[j, i] = this[j, i]
+        }
+    }
+}
+
 data class LMSettings (
     var iteration:Int,
     var func_calls: Int,

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

@@ -212,6 +212,36 @@ public fun DoubleTensorAlgebra.svd(
     return Triple(uTensor.transposed(), sTensor, vTensor.transposed())
 }
 
+public fun DoubleTensorAlgebra.svdGolubKahan(
+    structureND: StructureND<Double>,
+    iterations: Int = 30, epsilon: Double = 1e-10
+): Triple<DoubleTensor, DoubleTensor, DoubleTensor> {
+    val size = structureND.dimension
+    val commonShape = structureND.shape.slice(0 until size - 2)
+    val (n, m) = structureND.shape.slice(size - 2 until size)
+    val uTensor = zeros(commonShape + intArrayOf(n, m))
+    val sTensor = zeros(commonShape + intArrayOf(m))
+    val vTensor = zeros(commonShape + intArrayOf(m, m))
+
+    val matrices = structureND.asDoubleTensor().matrices
+    val uTensors = uTensor.matrices
+    val sTensorVectors = sTensor.vectors
+    val vTensors = vTensor.matrices
+
+    for (index in matrices.indices) {
+        val matrix = matrices[index]
+        val matrixSize = matrix.shape.linearSize
+        val curMatrix = DoubleTensor(
+            matrix.shape,
+            matrix.source.view(0, matrixSize).copy()
+        )
+        curMatrix.as2D().svdGolubKahanHelper(uTensors[index].as2D(), sTensorVectors[index], vTensors[index].as2D(),
+            iterations, epsilon)
+    }
+
+    return Triple(uTensor, sTensor, vTensor)
+}
+
 /**
  * Returns eigenvalues and eigenvectors of a real symmetric matrix input or a batch of real symmetric matrices,
  * represented by a pair `eigenvalues to eigenvectors`.