Reduced row echelon form - Gaussian Elimination

[epsi1on]

Hello, Is there any straight way to compute a "Reduced row echelon form" of a n*m matrix in CSparse.Net? Definition of Reduced row echelon form in wikipedia.

Thanks

[wo80]

Hello and sorry for delay!

CSparse cannot compute the RREF of a sparse matrix. Actually, I couldn't find a lot of info on this topic regarding sparse matrices. Do you need the RREF for solving or rank calculation? Have you tried QR?

[epsi1on]

Yes, i think so. What i need is to reduce the number of variables in a m*n eq. system where m is number of equations where a equations are linearly independently , and n is number of variables, from x0 to xn (unknown vector is named X). Now i can bind a variables into rest of variables.

Example:

[0 0 1 1 3 0 2]    [x0] = [3]
[0 0 2 6 1 0 5]  x [x1] = [1]
[0 0 3 7 4 0 7]    [x2] = [4]
                   [x3]
                   [x4]
                   [x5]

lets call it:

C . X = B

in this example m = 3, n = 7, a = 2 (because third equation is gained from first + second one, and is not a real equation that can reduce the number of variables), I want to select a=2 variables for example x1 and x2, and a matrix P, then write rest of variables regarding to these two.

Form a new matrix named CB is:

[0 0 1 1 3 0 2 | 3]  
[0 0 2 6 1 0 5 | 1]

RREF form of it (RREF of CB) is:

[+0.00 +0.00 +1.00 +0.00 +4.25 +0.00 +1.75 +4.25]
[+0.00 +0.00 +0.00 +1.00 -1.25 +0.00 +0.25 -1.25]

With this matrix, I can extract x2 and x3 from system:

[1 . x2 + 4.25 . x4 + 1.75 . x6 = 4.25]
[1 . x3 - 1.25 . x4 + 0.25 . x6 =-1.25]

Then after extract:

x2 = -4.25 . x4 - 1.75 . x6 - 4.25
x3 = -1.25 . x4 - 0.25 . x6 + 1.25

in matrix form:

[0 0 -4.25 0 -1.75]  [x0]   [-4.25]    [x1]
[0 0 -1.25 0 -0.25]  [x3]   [1.25]     [x2]
                     [x4] 
                   * [x5] +          = 
                     [x6]

Lets call it:

Cr * Xr = Br

where r means Reduced. In new system, variable count is reduced from n=7 to n-a=5. What I need to find is Cr matrix or CB matrix, i'm not sure if these can gain using any existing decomposition in the CSParse.NET

Thanks

[wo80]

CSparse.NET does not support this kind of reduction of the system. I took a closer look at SparseQR. It surely can solve your equations, but I'm afraid it won't help you finding the RREF.

What dimensions would your matrix usually have? You might want to write your own dense RREF code.

[epsi1on]

Is there a way to make my m*n matrix (where m < n) lower triangular?

maybe like this

[wo80]

Unfortunately, SparseLU doesn't support rectangular matrices and there are no plans to add more features to the sparse solvers of CSparse.NET (it's supposed to stay on par with the original CSparse).

I've been working a lot with native solvers lately, especially with SuiteSparse (Cholmod/Umfpack/SPQR etc.), and I'm writing C# bindings for them, which I plan to publish here on Github (maybe at the end of this month).

SPQR for example can easily extract the desired RREF (I just tried). Would the use of native DLLs be an option for you?

[epsi1on]

Thanks for reply, I rather to stay managed rather than native.' Another question: Suppose i have a m*n matrix with rank of o where m < n and o < m which means n-o rows are not independent and they are sum of multiply of other vectors. Is there any way to extract o independent row vectors of this system with current features of CSparse.net?

Thanks

[wo80]

You could try SparseQR and look at the R factor. Use reflection to get the non-public member:

CompressedColumnStorage<double> GetFactorR(SparseQR qr)
{
    var info = typeof(SparseQR).GetField("R", BindingFlags.Instance | BindingFlags.NonPublic);

    return info.GetValue(qr) as CompressedColumnStorage<double>;
}

Then, if abs(R[i, i]) > tol, you found an independent row. Make sure to take row permutaions into account (again, you can get the SymbolicFactorization containing row and column permutations using reflection):

SymbolicFactorization GetSymbolicFactorization(SparseQR qr)
{
    var info = typeof(SparseQR).GetField("S", BindingFlags.Instance | BindingFlags.NonPublic);

    return info.GetValue(qr) as SymbolicFactorization;
}

And be aware that SparseQR will factorize the transposed matrix since rows < columns.

[epsi1on]

thanks for detailed help. can you please also show then cpp code that computes RrEF form of matrix with SPQR?

thanks

[wo80]

I'm exporting the C API directly from SPQR, then do the interop call and the rest of the calculation in C#. It's way too much to post here. The function to call is named SuiteSparseQR_C. For your above matrix you'd get

Matrix R: 3x7
       0.0       0.0    -3.742    -9.087    -4.543       0.0    -8.820
       0.0       0.0       0.0    -1.852     2.315       0.0    -0.463
       0.0       0.0       0.0       0.0       0.0       0.0       0.0

Matrix Z (corresponding to rhs B): 3x1
    -4.543
     2.315
     0.000

Then find the pivot columns (dependent variables = left-most non-zero column (if any) in each row) (2 3 -1) and extract the corresponding submatrices of dependend and independent columns:

Matrix Rdep: 2x2
    -3.742    -9.087
       0.0    -1.852

Matrix Rindep: 2x5
       0.0       0.0    -4.543       0.0    -8.820
       0.0       0.0     2.315       0.0    -0.463

Solving the triangular system Rdep for Rindep then gives

Matrix X: 2x5
       0.0       0.0     4.250       0.0     1.750
       0.0       0.0    -1.250       0.0     0.250

which is the system you are looking for.

[epsi1on]

(again, you can get the SymbolicFactorization containing row and column permutations using reflection)

OK, the code to extract independent vectors is it? if qr.r[i,i] == 0, then row number qr.s.q[i] in original matrix is depended vector, otherwise row number qr.s.q[i] in original matrix is in-depended.

is this statement right?

[wo80]

Yes, this should be correct.

For a matrix A with rows < columns, SparseQR does a factorization of the transpose A', so S.q stores the column-ordering and S.pinv stores the row permutation of A'. Since you are interested in the original (non-transposed) system, the roles of the permutation vectors are reversed.

Take a look at the Solve method to see how the permutations are applied.

Though I don't know what exactly you are trying to do, I can't see how this will help you in reducing the system to anything like RREF. Even if you extract some independent rows, where to go from there?

And note: you probably don't mean qr.r[i,i] == 0. You will have to use a tolerance, which in turn means that you only get an approximation of "independent" rows.

[epsi1on]

thanks for information. Do you think it is easy or not for me to implement a class for finding RREF of matrix using Gaussian elimination operations regarding current state of CSparse and tools already in it? can you please show a free entry point for implementing such a class to me? (like pdf etc.) I need it so much at the time.

Thanks