A challenging problem: how to obtain a semi-lower triangular matrix from a general matrix?

Dear All,

I want to convert a general sparse matrix into a semi-lower triangular matrix. For example, a given matrix A:

A=[0 0 1 -1; -1 2 1 0; 0 2 -1 -1; -1 0 0 1;-1 -1 3 -1; 0 0 1 0; 0 0 -1 1;2 0 -1 -1];

Re-order the rows and columns, A can be converted into B:

B=[1 0 0 0; 1 -1 0 0; -1 1 0 0; 0 1 -1 0; -1 -1 2 0; 1 0 -1 2; -1 -1 0 2; 3 -1 -1 -1];

My code works well but it is very slow. For a 10000 by 10000 matrix, it takes more than 1 hour. For your information, I would like to share with you my algorithm.

My algorithm:

1. For a given sparse matrix A, we count the # of non-zeros elements in all rows, and re-order the rows from top to the bottom according to their # of non-zero elements (descend).

2. Deal with the first row whose # of non-zero elements is 1 (if a row with # = 1 does not exist, we select a row for its # = 2). If this "1" occurs at column j, switch the column 1 and column j.

3. Define the sub-matrix A(2:end,2:end) as a new A and repeat the above steps until the entire matrix is done.

The above steps will lead us to matrix B.

I do not know if someone could help find a faster algorithm than mine.

So far I got replies from Bruno Luong, OCDER, Matt J, and Stephen Cobeldick. Thank them all for their big help. But this equation is still open for solution. I hope you continue to investigate this challenging and interesting problem. I am looking forward to any better ideas.

Benson

Texas, USA

[Adam Danz]

tril() and it's cousin triu() will give you the trianulgar parts of a matrix. But I can't quickly see how you converted A to B. What were the conversion rules you used?

[Stephen23]

@Benson Gou: what is the exact algorithm to get from A to B?:

A =
   0   0   1  -1
  -1   2   1   0
   0   2  -1  -1
  -1   0   0   1
  -1  -1   3  -1
   0   0   1   0
B =
   1   0   0   0
   1  -1   0   0
   0   1  -1   0
   1   0  -1   2
  -1  -1   0   2
   3  -1  -1  -1

[Jos (10584)]

What is the rationale behind your re-ordering procedure?

[Bruno Luong]

I believe the problem is to find permutation of rows and columns of the original matrix to get as close as possible to a lower triangular matrix.

[Adam Danz]

The problem is not defined well enough to provide useful feedback. For example, are you just looking to take any of the non-zero values in A and just put them in any order filling the lower, left corner of B?

[Benson Gou]

Thanks you all for your responses. I would like to add more descriptions on my question: I want to convert a general sparse matrix (non-zeros elements are integers, and the sum of each row is zero except one row which only has one non-zero 1). Please verify it in matrix A given above. I want to add two requirements on the new matrix B: 1. Just like lower triangular matrix, each row only introduces one new non-zero column comparing with the previous rows. For our case, the difference is that we may have multiple rows introduce the same new non-zero column. For example, B=[1 0 0 0 0 0 1 -1 0 0 0 0 2 -1 -1 0 0 0 -1 0 1 0 0 0 0 0 1 -1 0 0 2 -1 0 -1 0 0 3 -1 0 0 -1 -1 0 0 -1 0 2 -1]; In matrix B, the 3rd and 4th rows introduce the same one non-zero column (column 3), and 5th and 6th rows introduce the same one non-zero column (column 4), and the last two rows introduce two new columns (columns 5 and 6). In another word, we do not want two new columns occur earlier than necessary (each row introduces only one non-zero column as possible as it can). We want to push them (2 new non-zero columns or more new non-zero columns) to the right of B. 2. Besides the above requirement, we hope we could push columns with most non-zero elements to the left of B. This requirement may be related to requirement 1.

If it is still not clear, please feel free to let me know. Thank you all for your great help with this interesting problem.

[OCDER]

Just curious, what is the application of this? And is it anything related to this : https://www.mathworks.com/help/matlab/ref/lu.html ?

If you have a working algorithm that's slow, upload that code. It's easier for us to rework it to be faster than to develop an algorithm from scrap that may not work as intended.

[Benson Gou]

Hi, OCDER,

My problem is not a LU decomposition. But thanks for your reply.

I tried your updated code, it is very fast and got a close answer. In the matrix B, your code gives: B(1:4,1:4) = [1 0 0 0;0 -1 1 0;0 1 0 -1;0 1 0 0];

What I hope the B matrix could be? B(1:2,1:4) = [1 0 0 0; X X 0 0].

In my updated description, I said I hope each new row could introduce only 1 new non-zero column (we keep this rule as possible as we can because sometimes we have to introduce 2 new non-zero columns, but we should keep it occurs closer to the right of matrix B). But your B, the second row introduces 2 new non-zero columns.

My code is very long. My idea is to re-order matrix A one row and one column at a time, so it is very slow. For a 10000 by 10000 it takes more than 1 hour.

Thanks a lot again for your great help.

[Bruno Luong]

I must admit I don't have courage to do through the description, you speak of 'push', 'introduce', 'new', etc..

It seems to me you have some sort of specific algorithm in your mind. If you could describe what is the desired shape of the final matrix B, it would be easier for us (me) to follow.

[Benson Gou]

Hi, Bruno, Thanks for your reply. I do have an algorithm to get B in my mind, but I think it is very slow.

My algorithm:

1. For a given sparse matrix A, we count the # of non-zeros elements in all rows, and re-order the rows from top to the bottom according to their # of non-zero elements (descend).

2. Deal with the first row whose # of non-zero elements is 1 (if a row with # = 1 does not exist, we select a row for its # = 2). If this "1" occurs at column j, switch the column 1 and column j.

3. Define the sub-matrix A(2:end,2:end) as a new A and repeat the above steps until the entire matrix is done.

The above steps will lead us to matrix B.

I do not know if you or someone else could help find a faster algorithm than mine.

Thanks a lot.

Benson

[Stephen23]

@Benson Gou: please do not close questions that already have answers.

[Benson Gou]

Hi, Stephen,

I opened this question again. But I improved my description. So far I did not get the right answer to my question. But I would like to thank you all for your great efforts.

Bei

[Stephen23]

@Benson Gou: please do not close questions that already have answers.

It actually makes it harder for us to help you, if you keep asking identical questions and closing them: it is harder for us to know what questions you have been asked, what information you have already given us, what answers you have been given, and what has been discussed already.

[Matt J]

@Benson,

It seems to me that your algorithm fails with the following input

A =
     0     0     0     0
     0     1     0     0
     0     0     1     1
     0     0     1     1

In step 3, the algorithm must eventually reduce A to its sub-matrix A(3:end,3:end)=ones(2) with no way to make that sub-matrix lower triangular.

However, it is clearly possible to make the original A lower triangular just by interchanging columns 1 and 4,

B =
       0     0     0     0
       0     1     0     0
       1     0     1     0
       1     0     1     0

[Benson Gou]

Dear OCDER and Bruno,

I run your code on the matrix A given as below: A=[0 0 1 -1; -1 2 1 0; 0 2 -1 -1; -1 0 0 1;-1 -1 3 -1; 0 0 1 0; 0 0 -1 1;2 0 -1 -1];

But your codes generated the following results:

It is close but still a little problem: 4th row and 7th row are not correct.

Thanks a lot for your efforts.

Bei

[Matt J]

@Benson,

Are you sure the thing you are ultimately trying to accomplish would not work just as well if you could triangularize T1*A*T2 for some pre-/post transformations T1, T2? You still haven't told us what you plan to do with the triangularized result.

Answer 1

[OCDER]

A = [ 0   0   1  -1;
     -1   2   1   0;
      0   2  -1  -1;
     -1   0   0   1;
     -1  -1   3  -1;
      0   0   1   0]; 
[~, SortedRowIdx] = sort(sum(A == 0, 2), 'descend');
B = A(SortedRowIdx, :);
[~, SortedColIdx] = sort(sum(B == 0, 1), 'ascend');
B = B(:, SortedColIdx);

 B =
     1     0     0     0
     1    -1     0     0
     0     1    -1     0
     1     0    -1     2
    -1    -1     0     2
     3    -1    -1    -1

Is the issue that this is too slow though?

[Matt J]

I don't think it's reliable enough. With this matrix, it fails.

A = [ 0   0  -1   0;
     -1   2   1   0;
      0   2  -1  -1;
     -1   0   0   1;
     -1  -1   3  -1;
      0   0   1   0];

The proposed algorithm gives

   B =

      -1     0     0     0
       1     0     0     0
       0    -1     0     1
       1    -1     2     0
      -1     0     2    -1
       3    -1    -1    -1

This is not upper triangular, but will be if you do one more row interchange (rows 3 and 4).

[Bruno Luong]

The question is can you do better?

[Matt J]

Yes, you can do better. Just interchange rows 3 and 4.

[Bruno Luong]

Sorry I don't see it. interchange rows of what? B "incorrect"?

And also what criteria tells "B correct" is more correct than the other? Both have a cluster of 7 zeros.

[Matt J]

Yes. Run the proposed algorithm to obtain,

    B =

      -1     0     0     0
       1     0     0     0
       0    -1     0     1
       1    -1     2     0
      -1     0     2    -1
       3    -1    -1    -1

Now just interchange row 3 and 4 of this matrix and it will become upper triangular.

[Benson Gou]

Dear Bruno,

Thanks for your great help. Your code is very fast. But there is one thing I want to point out: like lower triangular matrix, each row introduces a new non-zero column, I hope each new row only introduces a new non-zero column comparing with previous rows. Please see the new description of my problem.

Thanks a lot again. Benson.

[Benson Gou]

Dear OCDER,

Thanks lot for your great help. I run your code but did not get a right matrix B. Would you please check your code again?

Thanks a lot again. Benson