Translation from NumPy to EJML

[CLOVIS]

Hi,

I am new to EJML, and I'm trying to rewrite a file from NumPy.

I'm currently lost on the meaning of the following line:

result, residual, rank, s = numpy.linalg.lstsq(self.A, self.b)

where A is a 3*n matrix (which, if I understood correctly, should be a FMatrixRMaj), and b is a vector of size 3.

I've been reading the 'Solving Linear Systems' documentation page and it mentions the LinearSolverFactory_DDRM.leastSquares() function, which sounds like it's what I'm searching for, but I have no idea how to actually use it.

For additional context, the file I'm trying to rewrite is this one.

Thanks for any help ^^

[CLOVIS]

Here is my current progress:

After analyzing the original code, I noticed that 'rank' and 's' are unused, so that's already easier.

I managed to get the solution using:

var solver = LinearSolverFactory.FDRM.leastSquares(A.numRows, A.numCols);

solver.setA(A);

var solution = new FMatrixRMaj(A.numRows, A.numCols);

solver.solve(b, solution);

According to my unit tests, that is in fact the results I'm expecting. I just lack the 'residual' part for now.

[Peter A]

You might want to consider starting off using `CommonOps_DDRM.solve()` just because the API is less verbose. Out of curiosity why did you choose to go with the low level API on that page? You probably also want to use 64-bit floats, which means using DMatrixRMaj and stuff with the suffix `DDRM`, but maybe you're on a very low memory system.

Looks like numpy is doing a bit of extra work here that's not standard. There's some abuse of definitions on NumPy's side but I won't go into that. They define the residual as `|b-A*x|` in NumPy's manual. So what you need to do is multiply `bp=A*x` then compute `r=b-bp`, and finally call `|b-A*x| = NormOps.normf(r)` to get your solution. You can avoid the creation of one of the temporary matrices by using this function that's in a more obscure corner of the code `|b-A*x| = SpecializedOps_DDRM.diffNormF_fast(b,bp)`.

--
You received this message because you are subscribed to the Google Groups "efficient-java-matrix-library-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to efficient-java-matrix-li...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/efficient-java-matrix-library-discuss/b32b6e8e-b8d7-4c21-9db9-47bc36f12ebbn%40googlegroups.com.

--

"Now, now my good man, this is no time for making enemies."    — Voltaire (1694-1778), on his deathbed in response to a priest asking that he renounce Satan.

[CLOVIS]

Thanks a lot, I'll try to write it with this API :) Honestly it says the low level API is a bit better when it comes to performance, and since I'm expecting this bit of code to be ran very often, I'm trying to put all the chances on my side—and I do realize it might be a bit naive to try and do that when I know close to nothing about the library and I don't have a real benchmark to back this up, but well...

For this use case I don't really need the extra precision of 64-bit floats so I figured I'd get a little out of it, but speed is more important than memory usage here, so maybe you have some information on what to choose?

Thanks for the given functions, I managed to write something that works, but it's probably not that performant, I'll rewrite it with your functions:

// Solving A x = b
var solver: LinearSolverDense<FMatrixRMaj> =
LinearSolverFactory_FDRM.leastSquares(normals.numRows, normals.numCols)
solver.setA(normals)
val solution = FMatrixRMaj(normals.numRows, normals.numCols)
solver.solve(b, solution)

// Calculates b - A x
var rs2 = new FMatrixRMaj(normals.numRows, normals.numCols)
MatrixVectorMult_FDRM.mult(normals, solution, rs2)
SimpleOperations_FDRM().minus(b, rs2, rs2)

// Sums the square of the results, to get the residual
var residual = 0f
for (int i = 0; i < rs2.getNumElements(); i++)
    residual += rs2.get(i) * rs2.get(i)