Linear system Ax=b solver for Matrix abstraction

[LoicBer]

When working with the matrix abstraction, the user may have to solve Ax=b linear systems where A is a square matrix and b a vector. For example, in the Alternating Least Squares algorithm with k features, each step requires multiple Ax=b resolutions where A is a k*k matrix. Since k is not supposed to be a large number, these calculations could be done on a single node using direct solving methods based on matrix factorization (Cholesky LL, LDL, LU, PLU...).

So, I implemented a %\% Matrix operator where A %\% b computes x, the solution of Ax=b without computing explicitly the inverse of A. The proper method to solve Ax=b depends on the properties of A. For example, if A is already triangular, there is no need to run a heavy factorization method like PLU. Thus, the solver use the following work flow.

Notes :

• It is inspired by the Matlab work flow for Ax=b solving which you can find here http://fr.mathworks.com/help/matlab/ref/mldivide.html. Simplified because we do not have complex matrices.
• These are only direct solving methods. They are full of dirty dependencies which make them hard to parallelize, they suit well for matrices that can fit in a single node.
• We could also imagine to apply these methods to the blocks of a larger matrix and then combine the results. For example the inverse of a matrix can be done by block. But I do not know yet if Ax=b can be solved by block without explicitly computing inv(A).
• If you wonder why the Cholesky method may fail (by leading to the computation of a square root of a negative number), it is because it is only guaranteed to work on symmetric positive-define matrices. But the positive-define criterion is really hard and expensive to check.

For the moment, the 3 factorization methods (in gray), the triangular solvers (in yellow) and the whole Ax=b solver are public methods of the Matrix class. For instance
val m: Matrix[Double] = ...
val (p,l,u) = m.pluFactorization()
returns the 3 matrix factors of the PLU decomposition of m.

Do you think that we should provide these factorization methods to the user or should we just make them all private and only provide the whole %\% solver (because no one would ever need to get the matrix factors for an other task ) ? This also raises the question of Matrix properties tests. For the moment, the properties are checked by the general solver in order to choose the right factorization method. But, If the factorization methods are public, should they also check that the Matrix have the right properties ? For now, they have also some assert(...) to check that the computation makes sense. But It means running twice the same tests in the general solver. Should we add an optional boolean parameter to run or not the tests in the factorization methods ?

What could also be interesting would be to have an object MatrixUtils with methods to check matrix properties, like isUpperTriangular(A: Double): Boolean, and more complex properties... We could also put the factorization methods there instead. What do you think ?