#clapack_gesv and #cblas_trsm for matrix-valued systems require taking transposes before and after

[agisga]

For my linear mixed models gem I was trying to figure out a way to solve linear systems of the form A*X = B, where A is a nxn matrix, B and X are nxp matrices (without taking the inverse of A). After much trying, the only way that works is

bt = b.transpose
NMatrix::LAPACK::clapack_gesv(:row, a.shape[0], b.shape[1], a, a.shape[0], bt, b.shape[0])
x = bt.transpose

That is, I have to take an additional transpose of b before applying clapack_gesv, and then transpose the result.
Interestingly, that issue does not happen when b is a vector (i.e. no extra transposing required in that case).

Here is an example calculation that I tried out (an example of a linear system with one-dimensional b that does not require any #transpose can be found in the spec):

a = NMatrix.new([3,3], [1, 2, 3, 0, 0.5, 4, 3, 3, 9], dtype: :float64)
b = NMatrix.new([3,2], [1, 1, 2, 2, 3, 3], dtype: :float64)
bt = b.transpose
NMatrix::LAPACK::clapack_gesv(:row,a.shape[0],b.shape[1],a,a.shape[0],bt,b.shape[0])
x = bt.transpose
a = NMatrix.new([3,3], [1, 2, 3, 0, 0.5, 4, 3, 3, 9], dtype: :float64)
a.dot x
# => 
[
  [0.9999999999999999, 0.9999999999999999]
  [               2.0,                2.0]
  [2.9999999999999996, 2.9999999999999996]
]

It does not work without taking transposes. Or am I doing something else wrong?

[agisga]

I get the same issue when trying to solve triangular matrix-valued systems with #cblas_trsm.
If a is lower triangular, and b is a nx1 NMatrix the following solves the system A*x=b just fine:

NMatrix::BLAS::cblas_trsm(:row, :right, :lower, :transpose, :nounit, 
                              b.shape[1], b.shape[0], 1.0, a, a.shape[0],
                              b, b.shape[0])

But if b is a nxm NMatrix with m>1, then I have to take transposes:

bt = b.transpose
NMatrix::BLAS::cblas_trsm(:row, :right, :lower, :transpose, :nounit, 
                              b.shape[1], b.shape[0], 1.0, a, a.shape[0],
                              bt, b.shape[0])
x = bt.transpose

I cannot make it work without taking transposes...

[wlevine]

Argh, I just went through and tested all these functions. I actually thought gesv was working fine, but only because I only tested the case where b is a vector.

[translunar]

You two are running into some of the problems I encountered when testing these functions. There are a lot of functions, each with edge cases. My apologies for creating bugs years ago that are now haunting you. =/

[atimin]

I'll try to fix it.

[atimin]

I've just got around to the bug. Above all NMatrix::BLAS::cblas_trsm is ok. @agisga, you called the method with wrong arguments. This is right for A*X = B:

NMatrix::BLAS::cblas_trsm(:row, :left, :lower, :transpose, :nounit, 
                              b.shape[0], b.shape[1], 1.0, a, a.shape[0],
                              bt, b.shape[0])

I compared it with C-function from cblas. Their behaviours are identical.

I had a look at NMatrix::LAPACK::clapack_gesv too. Before I start working on it I'd like to know why we don't use the ATLAS's function? It might be better to rewrite the method at all.

[agisga]

That doesn't seem to work for me (assuming the bt in your post is a typo and should be a b):

[52] pry(main)> a
=> 
[
  [1.0, 0.0, 0.0]   [2.0, 3.0, 0.0]   [4.0, 5.0, 6.0] ]
[53] pry(main)> b = NMatrix.new([3,3],(1..9).to_a, dtype: :float64)
=> 
[
  [1.0, 2.0, 3.0]   [4.0, 5.0, 6.0]   [7.0, 8.0, 9.0] ]
[54] pry(main)> NMatrix::BLAS::cblas_trsm(:row, :left, :lower, :transpose, :nounit, 
[54] pry(main)*   b.shape[0], b.shape[1], 1.0, a, a.shape[0],                              
[54] pry(main)* b, b.shape[0])  
=> true
[55] pry(main)> a.dot b
=> 
[
  [-2.4444444444444446, -2.2222222222222223, -2.0]   [ -6.722222222222223,  -6.111111111111111, -5.5]   [ -5.833333333333336,  -3.666666666666666, -1.5] ]
# This should be [ [1.0, 2.0, 3.0]   [4.0, 5.0, 6.0]   [7.0, 8.0, 9.0] ]

Compared to my workaround from above, which works fine:

[56] pry(main)> a
=> 
[
  [1.0, 0.0, 0.0]   [2.0, 3.0, 0.0]   [4.0, 5.0, 6.0] ]
[57] pry(main)> b = NMatrix.new([3,3],(1..9).to_a, dtype: :float64)
=> 
[
  [1.0, 2.0, 3.0]   [4.0, 5.0, 6.0]   [7.0, 8.0, 9.0] ]
[58] pry(main)> bt = b.transpose
=> 
[
  [1.0, 4.0, 7.0]   [2.0, 5.0, 8.0]   [3.0, 6.0, 9.0] ]
[59] pry(main)> NMatrix::BLAS::cblas_trsm(:row, :right, :lower, :transpose, :nounit, 
[59] pry(main)*   b.shape[1], b.shape[0], 1.0, a, a.shape[0],                              
[59] pry(main)* bt, b.shape[0])                                
=> true
[60] pry(main)> x = bt.transpose
=> 
[
  [                  1.0,                  2.0,  3.0]   [   0.6666666666666666,   0.3333333333333333,  0.0]   [-0.055555555555555504, -0.27777777777777773, -0.5] ]
[61] pry(main)> a.dot x
=> 
[
  [1.0, 2.0, 3.0]   [4.0, 5.0, 6.0]   [7.0, 8.0, 9.0] ]

[atimin]

Sorry, I wrote the wrong call, it's necessary to use false instead of :transpose in line 54.

NMatrix::BLAS::cblas_trsm(:row, :left, :lower, false, :nounit, 
                              b.shape[0], b.shape[1], 1.0, a, a.shape[0],
                              b, b.shape[0])

the method is really unfriendly because it's hard to write a call without errors.

[agisga]

Thanks!
Also the last argument in the above should be b.shape[1] (otherwise it only works on square matrices b).
Should the spec be changed to display the correct function call?

[atimin]

I think, yes. I'll do it.

But what about NMatrix::LAPACK::clapack_gesv? My answer is still unanswered. @MohawkJohn could you make it clear?

[translunar]

Honestly, I can't recall. It's been too long. I thought it was using the ATLAS function.

[wlevine]

My thought is that we should just get rid of clapack_gesv and just make #solve work for non-vector equations.

[wlevine]

It seems like when you call ATLAS clapack_dgesv directly from C, you get the same behavior. I think what is happening is that when you specify CblasRowMajor, it only ends up interpreting the first argument a as row-major, but still treats b as if it were column-major. So I think for nmatrix transposing b before and after is really the best workaround we can come up with.

[wlevine]

Just for reference, here's the code I used to test clapack_dgesv: https://gist.github.com/wlevine/e0ee4a9a0f527b666c65

[wlevine]

Ah, here's the entry in the ATLAS FAQ explaining this behavior: http://math-atlas.sourceforge.net/faq.html#RowSolve

[agisga]

@wlevine thanks for figuring it out. In my opinion the best solution in this case is your idea of making #solve work for linear systems with a matrix as right hand side.

Regarding the unfriendliness of trsm, I can open a pull request with the #triangular_solve method, that I use in mixed_models. Should I do that?

[wlevine]

@agisga In the current version of nmatrix, I've removed clapack_gesv and made #solve work with non-vector right-hand sides. If you can verify that #solve works for you, you can close this bug.

I think a pull request with #triangular_solve would definitely be welcome.