RFC: factorize() and \

[andreasnoack]

A new general matrix factorization method. Some changes to the treatment of symmetric and hermitian matrices and some bug fixes. And some tests.

This is a follow up on the discussion in JuliaLang/LinearAlgebra.jl#16. The main new thing is factorize which tries to guess a good factorization depending on the structure of matrix. For example

julia> A=randn(5,5);

julia> typeof(factorize(A))
LU{Float64}

julia> typeof(factorize(A + A'))
BunchKaufman{Float64}

julia> typeof(factorize(A'A))
Cholesky{Float64}

julia> typeof(factorize(A[:,1:2]))
QRPivoted{Float64}

For complex indifinite matrices it is necessary to distinguish between symmetric and hermitian matrices (the Bunch-Kaufman case) and therefore I have added a Symmetric type which also works for real matrices. Hence Hermitian is now a complex matrix only type. I have also added the bkfact method for constructing the Bunch-Kaufman decoposition directly.

To handle least squares, I had to reimplement the functionality of xgelsy in Julia because I wanted to be able to solve rank deficient problems. You can now reuse a QR based "predicter":

julia> mX = randn(10,2)*randn(2,3);y1 = randn(10);y2 = randn(10);

julia> xf = factorize(mX);

julia> xf\y1
3-element Float64 Array:
 -0.18619   
  0.811663  
  0.00219657

julia> xf\y2
3-element Float64 Array:
  0.0448524
 -0.197673 
  0.0139636

I got some random failures in the tests but I couldn't get them if I ran the tests with include from a Julia session.

[StefanKarpinski]

This is really cool! I'm gonna let some of the numerical linear algebra experts take a look though.

[stevengj]

It seems more sensible to have both Hermitian and Symmetric work for both real and complex matrices, with Hermitian(A) equivalent to A+A' and Symmetric(A) equivalent to A+A.'. (So that Hermitian and Symmetric are equivalent for real matrices but inequivalent for complex matrices.)

[stevengj]

I remember that a colleague of mine once checked, and LU was actually faster than Bunch-Kaufman for symmetric-indefinite matrices, although the latter requires less memory. Is this still true?

[andreasnoack]

@stevengj Thank you for the comments. I think you are right about the Symmetric/Hermitian thing. I hadn't tested LU vs Bunch-Kaufman before your post but expected the LAPACK defaults to be good choices. Based on a few tests, it seems that for the complete solution of an indefinite problem, i.e. factorize+solve, Bunch-Kaufman is fastest. However, the factorization is not as clean as the LU and therefore each application of the factorization will be slightly slower. I haven't looked into the relative memory consumption of the two procedures.

[ViralBShah]

We should perhaps also check for Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal. In the triangular systems, we also need to do checks to look for permutations of triangular matrices.

I wonder if we can rewrite things to classify a matrix in the minimal number of sweeps over the matrix.

[ViralBShah]

Let's merge this as soon as it is stable, and more intelligence for various matrix types can be added over time.

[andreasnoack]

@ViralBShah I tried to write a structurize function which checked for all the types you mentioned. It is very likely that it can be rewritten to become faster but I fear that too much time will be spend on checking the structure relative to the time gain of using a specialized solver.

@dmbates It would be good if you could have a look at the \(QRPivoted,Matrix) method at line 306 in factorization.jl. It is very central to the least squares functionality in Julia after this pull request is merged. Also, now that the rank detection is done in Julia, maybe it should be part of the pivoted QR factorization. In the present implementation, the rank would have to be determined for each application of the QRPivoted factorization.

[ViralBShah]

The time may not matter for factorize, which gets called once and reused often. For single use like \, in the case of small problems, this could be an issue. I say let's get it all in, and we can optimize later.

[andreasnoack]

That is right. I'll finish the changes now such that we can get this merged.

[dmbates]

@andreasnoackjensen The Julia implementation of xgelsy looks very good. I appreciate your taking the initiative.

A couple of things occur to me. I'm not sure of the motivation of your check on abs(r[1}) < rcond. The condition number of a matrix is a relative criterion. That is, if you multliplied the original model matrix by 1e-50 you shouldn't change the condition number but you could make it fail this test. Are you sure that test makes sense?

Two parenthetical remarks here - your code calculates abs(r[1]) three different times; you should only do that once (the compiler may recognize this but it is better to avoid this in the code). Secondly I saw an idiom in code by Kevin or Jameson of writing checks on argument values like

abs(r[1]) >= rcond || error("Left hand side is zero")

which I rather like. I often use the stopifnot function in R for this purpose.

And while I am being pedantic and picky, the error message in tzrzt! from lapack.jl should say "fewer columns".

One thing about calling xgelsy is that the transformed response Q'y is available upon completion. Occasionally it is an advantage to have that because the residual sum of squares can be calculated from that without evaluating the fitted values or the residuals. In fact, the sequential analysis of variance terms can be evaluated from that too.

Is it necessary to extract r, which will involve creating a copy? Most of the time LAPACK code can work with the full rectangular matrix as if it were square and upper triangular? It may make the code uglier but if this is a core computational routine it can be worthwhile avoiding copies.

[andreasnoack]

@dmbates Thank for the comments.They are appreciated. I'll follow up on them as soon as possible.

[dmbates]

@andreasnoackjensen You're welcome. And my thanks to you for all the work you have put into the linear algebra facilities. I think the community and especially the core developers are building something quite marvelous here!

[ViralBShah]

Some parallel test failed for clang, but gcc passed.

[andreasnoack]

I think this one is about to be ready to get merged. The suggestions of @dmbates have been incorporated except for the option of returning Q'y. I think that part should be handled when solve! methods are written as proposed in #3239. In the present implementation \ does not call factorize because the checking in factorize is only worth the time if the factorization is reused. However, \ still checks if the matrix is triangular. I have disabled a test because of #3396.

[andreasnoack]

I haven't touched the parallel stuff so I think it must be a Travis thing.

[ViralBShah]

Please merge when you are ready to do so. I think it is ok to not have \ call factorize right now, but along with the triangular systems, I would perhaps want to include a diagonal matrix check in \. Eventually, if a matrix is large enough, it will be faster to call factorize in \, but I think that can be the subject of a future PR.

[fredrikekre]

Was there a reason for the fallbacks on this and the next line, allocating a new array rather than calling BLAS.symm!?
Only the fallback from the line below are still on master, should that fallback be removed? How would you else be able to dispatch here or here
cc @andreasnoack

[andreasnoack]

I'm not completely sure but I think OpenBLAS' symmetric kernels used to be quite a bit slower than gemm and maybe that is the reason for this. We should revisit and probably use symm if it is not slower.