LU factorization returns both lower and upper triangular matrices in the same matrix

[CaptainSharf]

Summary of Problem

LU factorization returns both lower and upper triangular matrices in the same matrix.

Steps to Reproduce

Source Code:

use LinearAlgebra;
var dom = {1..3,1..3};
var arr:[dom] real;
arr[1,1] = 1.0;
arr[1,2] = 2.0;
arr[1,3] = 3.0;
arr[2,1] = 4.0;
arr[2,2] = 5.0;
arr[2,3] = 6.0;
arr[3,1] = 8.0;
arr[3,2] = 10.0;
arr[3,3] = 12.0;

var x = lu(arr);
writeln(x[0]);

Compile command: chpl LUTest.chpl

Execution command: ./LUTest

Output

8.0 10.0 12.0
0.5 0.0 0.0
0.125 inf nan

Configuration Information

• Output of chpl --version: 1.24.0 pre-release

[CaptainSharf]

@ben-albrecht , Based on the return type for this procedure in other popular frameworks, would it be a good idea to return L,U and P as separate entities?

[ben-albrecht]

Based on the return type for this procedure in other popular frameworks,

Do you have some specific examples?

I believe returning a single array was a memory-optimization in the UI. We can consider changing it if there is enough justification. If we left it as is, I think we could add a note to the documentation to show how one would get the lower/upper part of the matrix via tril and triu.

[bradcray]

I'm not a L.A. expert, but I've been under the impression that some scientific libraries store L and U in a conjoined fashion like it seems we do. We could potentially support both modes through a param argument or variant routine.

[ben-albrecht]

@bradcray - yes, I think libraries do both. As 1 data point, I see that scipy.linalg.lu returns (p, l, u) separately.

We could potentially support both modes through a param argument or variant routine.

That could be useful, though I generally like to avoid args that change the return type unless there's strong justification.

Just to write it out for reference, here's what unpacking the lower/upper result looks like today:

var (A, p) = lu(arr);
var lower = tril(A),
    upper = triu(A);

Addendum: This is for square matrices only, as @CaptainSharf points out below.

[CaptainSharf]

Hi @ben-albrecht , Thanks for giving the example! I too had scipy.linalg.lu in mind. In general, when matrices are rectangular, we can get upper using triu(A), however to get the lower part, I think we need to do

var (m,n) = A.shape;
var Lrect = A-U+eye(A.domain);
lower = Lrect[1..m,1..m]

The native implementation of lu doesn't support rectangular matrices. I have made a few changes in this procedure locally to accommodate both square and rectangular matrices. I would like to raise a PR to validate my changes

[ben-albrecht]

The native implementation of lu doesn't support rectangular matrices.

I think that's worth filing a separate issue -- would you be interested in opening it?

In general, when matrices are rectangular, we can get upper using triu(A), however to get the lower part, I think we need to do

I see. Do you mind sharing a full stand-alone example that I can try out? I'm not sure I follow your snippet entirely.

[CaptainSharf]

Hi @ben-albrecht , I have added a few tests in the file testLU.chpl as part of #16980 to validate our native implementation of lu factoring. I have used l,u to reconstruct the original matrix. I think maybe that'll illustrate what I have in mind. Please do let me know if it's still unclear. Note that the test file uses blas for the dot procedure, so make sure to add -lblas flag to compile it.

[vrockz747]

Hey, @CaptainSharf, I would like to help on this issue, could you please explain where I should focus and what can I help with.

[CaptainSharf]

Hi @vrockz747 , I have created #16979 , #16980 for fixing this issue, you could take a look at that PR