Hello,

New to Eigen. The linear algebra problem I am trying to solve computationally is as follows:

Let A be an upper unitriangular matrix (ie upper triangular with only ones on the main diagonal) and entries in the field of rational numbers. The matrix A has only occasional non-zero entries above the diagonal and may be considered to be sparse.

Let B be an invertible matrix with at most two non-zero entries in each row and column, and entries also in the field of rational numbers.

The problem is to solve for the matrix X which satisfies A * X = B * X A * X = B * A (matrix products). This problem needs to be solved exactly, ie float or double approximations for the rational entries will not suffice. In fact, the matrix B is part of a finite sequence {B_i} and the problem requires computation of an X_i for each B_i in the sequence, with A being fixed throughout.

The matrices A and B can rapidly get large. While I've been able to compute examples with dimensions up to 3000 * 3000 using closed source mathematical software like Mathematica and Maple, these calculations have quickly become impracticably time consuming at dimensions above 2500.

So far, using the Boost rational library, I've been able to construct matrices of the form Eigen::Matrix<boost::rational<int>, Eigen::Dynamic, Eigen::Dynamic>. Matrix multiplication works fine with these objects.

My questions are:

1) How does one go about solving for X in the equations A * X = B * X A * X = B * A in Eigen given that the entries of A and B are not floats or doubles? Since A is already upper triangular, this really comes down to whether Eigen has a suitable method of forward substitution that can be applied.

2) Would there be any computational advantage to constricting A and B as sparse rather than dense matrices?

Thanks in advance.

Antonio

Edit: typos.

Just to be specific, it looks like A, B, and X are all n * n matrices, with A and B invertible.

2) A is described to be sparse. Using a sparse matrix type has the benefits of only storing and operating with nonzero entries (plus a little overhead). You should use a sparse type to store A. B is dense, so the memory overhead of making it sparse incurs a cost. So B should be a dense type. Sparse and Dense matrices are compatible in a lot of operations, so don't worry too much about mixing the types.

And are you sure the system is A*X = B*X? Because then you would be solving (B - A)*X = 0 (with 0 being the n by n zero matrix), .

1) Solving the (B - A)*X = 0 problem is equivalent to solving an n^2 by n^2 system (in the usual form of M*y = 0, with M being n^2 by n^2 and y and 0 being n^2 column vectors). Can you see how symbolically multiplying the first column of X, x1, by the first row of (B - A), r1, gives r1*x1 = 0, which defines the coefficients of the first row in M for M*y = 0? In general, a 3000 * 3000 system like A*X = B*X this is actually a 9000000 * 9000000 system (ouch). In addition, you need exact solutions. So the complexity of directly solving a massive system is not in the favor of your time.

And are you sure the system is A*X = B*X? Because then you would be solving (B - A)*X = 0 (with 0 being the n by n zero matrix), .

Thank you for your reply and apologies for the error. Can't believe I let that error slip through, not once but twice.

The system I meant to write down was A * X = B * A or equivalently X = A^{-1} * B * A. It really is a change of basis problem, with both A and B invertible.

Although the dimensions of the matrices are quite large, the fact that A is upper triangular with ones on the diagonal, with only occasional non zero entries above the diagonal, should make it more tractable than appears at first sight. I've successfully run these calculations in interpreted languages, on cruddy hardware, in a few hours.

Thank you for the suggestion re storing A as sparse and B as a dense type.

Antonio