[EJML] Solving matrices returns NaN's as a result

[Melodeiro]

Im trying to run the following code:

import koma.*
import koma.extensions.get
import koma.matrix.Matrix

fun main(args: Array<String>) {
    var index = 1.0
    val a = fill(3,3) { _, _ -> index++ }
    val b = zeros(3,1)
    val c = fill(3,1) { row, _ -> row.toDouble() + 1 }
    printMatrix(a)
    printMatrix(b)
    printMatrix(c)
    val x1 = a.solve(b)
    val x2 = a.solve(c)
    printMatrix(x1)
    printMatrix(x2)
}

fun printMatrix(matr: Matrix<Double>) {
    for (i in 0 until matr.numRows()) {
        for (j in 0 until matr.numCols()) {
            print("${matr[i, j]} ")
        }
        println()
    }
    println()
}

The problem is that it returns 2 matrices with 3 NaN's. Should be zeros for x1, and random values for x2. If A * X = B where B = 0, there should be at least 1 answer: X = 0. Matrix with zeros in that case.

Happens, when im using koma-core-ejml while koma-core-mtj working as expected. Pure EJML code works fine too.

[kyonifer]

Thanks for the report. I can reproduce the issue locally. Can you provide the pure EJML code that ran successfully for you? I'm seeing the following using pure EJML and the high-level solve interface:

    val a = SimpleMatrix(3, 3, true, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0)
    val b = SimpleMatrix(3, 1, true, 0.0, 0.0, 0.0)
    val out = SimpleMatrix(3,1, true, 0.0, 0.0, 0.0)

    a.solve(b) // SingularMatrixException

    CommonOps_DDRM.solve(a.ddrm,b.ddrm,out.ddrm)
    println(out) // NaN, NaN, NaN

So it seems that EJMLs solve cannot handle singular A matrices (just like numpy's solve cannot). Both EJML and koma provide a pseudoinverse function that can be used to solve for singular matrices. In pure EJML:

    val x = a.pseudoInverse().times(b)
    println(allclose(EJMLMatrix(x), zeros(3,1))) // true

Similarly for your other case where b=[1,2,3]:

    val b2 = SimpleMatrix(3, 1, true, 1.0, 2.0, 3.0)
    val x2 = a.pseudoInverse().times(b2)
    print(allclose(EJMLMatrix(a*x2), mat[1,2,3].T)) // true

Doing the same thing in koma with EJML backend using b=[1,2,3]:

    val aKoma = mat[1,2,3 end 4,5,6 end 7,8,9]
    val bKoma = mat[1,2,3].T

    val xKoma = aKoma.pinv() * bKoma
    print(allclose(aKoma*xKoma, mat[1,2,3].T)) // true

[Melodeiro]

Equation version:

    val eq = Equation()
    eq.process("A = [1 2 3; 4 5 6; 7 8 9]")
    eq.process("B = [0; 0; 0]")
    eq.process("X = B / A")
    val x = eq.lookupMatrix("X")
    println(x)

Oh, and also i just realised, that i used for SimpleMatrix version the same code: eq.process("X = B / A") and eq.alias() to put the matrices into equation. Maybe its an issue with EJML algo then.

[kyonifer]

eq.process("X = B / A")

Curious, according to their docs this should also throw an exception for singular matrices. However, I get the same result as you, it works.

At any rate, I do agree that the user experience would be better if the Matrix facade was consistent across backends, and also if it "just worked" on solving singular matrices.

One solution is for us to try calling EJMLs CommonOps_DDRM.solve and if we get a failure status on return we try using the pseudoinverse. It'd probably be more efficient if EJML has a codepath that doesn't involve trying to solve twice-- they must have some codepath solving via SVD if the equation version can solve singular matrices.

Another solution is for us to always call the pseudo-inverse. There can be differences with numerical conditioning between LU/QR and SVD solvers and typically SVD solvers are more expensive, though I haven't looked into benching EJMLs methods against each other in quite some time to know how much effect this will have.

The first solution is a quick fix for now anyway.