Solving Linear equation using LU Factorization and COLAMD

I am trying to solve the Poisson equation on a rectangular/square domain with homogeneous Neuman boundary conditions. When discretized using central finite differences, I end up with an equation Ax = b, where A is sparse and singular. A can be rendered invertible by the following modification A(1,:)=0; A(1,1)=fixed value (e.g., 1000), such that we specify the value of x at a location in the domain.Typically, the equation is efficiently solved by prefactoring A using LU factorization. i.e., x = U\(L\b).

By permuting the columns of A (using COLAMD), however, the fill-ins of the resulting L,U matrices can be decreased, which can speed up solving the equation. So, p = COLAMD(A), [L,U] = lu(A(:,p));x(p) = U\(L\b(p)). Here, p indicates the column permutations.

Here is the problem I have.Let x1 and x2 be the solutions computed without and with column permutations. I compute ||LU-A(:,p)|| <1E-10, to test the COLAMD is doing what it is supposed to do. However, the solutions x1 and x2 are very different. ||x1-x2||. I am pasting the code below. Please let me know.

nx = 10; ny = 10; 
% Construct matrix A
d2Qdx2 = spdiags(ones(nx,1)*[1, -2, 1],[-1,0,1],ny,nx);
d2Qdx2(1,1) = -1; d2Qdx2(nx,nx) = -1;
d2Qdx2 = kron(d2Qdx2,speye(ny));
d2Qdy2 = spdiags(ones(ny,1)*[1, -2, 1],[-1,0,1],ny,ny);
d2Qdy2(1,1) = -1; d2Qdy2(ny,ny) = -1;
d2Qdy2 = kron(speye(nx),d2Qdy2);
A = d2Qdx2 + d2Qdy2;
A(1,:)=0; A(1,1) = 1;

% generate the right hand side vector
B = rand(ny,nx); b = B(:)-mean(B(:)); 

% Solution without COLAMD
[L,U] = lu(A); x1 = U\(L\b);
% solution with COLAMD
p = colamd(A); [L,U] = lu(A(:,p)); x2(p) = U\(L\b(p));
norm(L*U-A(:,p),'fro')

ans = 3.5425e-15

norm(x1-x2,'fro')

ans = 74.5403

Answer 1

[John D'Errico]

First, making a 100x100 matrix sparse is silly. Flat out, silly. The amount of time needed to solve that system is tiny, even as a full matrix.

Is A no longer singular as you have formed it?

cond(full(A))
ans =
       1043.6

Not even close to singular. And that means any column permutations will have no impact on the solution.

timeit(@() A\B)
ans =
   7.2113e-05

It took fractions of a millisecond to solve. In full form, is it any different?

Af = full(A);
timeit(@() Af\B)
ans =
    7.734e-05

Virtually, NO. So unless your real problem is wildly larger, then you are wasting your time with this entire effort to make A sparse, or for that matter, to use column permutations to reduce fill-in. Even if you made this a 1000x1000 problem, the direct full solve will still be blazingly fast.

It would help if you tell us how large your real problem is, or why you so desperately want to solve this problem as you are doing it.

Even without that, you still clearly have a problem in how you employed the column permutaiton.

b = B(:)-mean(B(:)); 
[L,U] = lu(A); x1 = U\(L\b);
norm(x1 - A\b)
ans =
     0

x1 is computed correctly. Now let me look at how you are computing x2, using the column permutation.

What does a column permutation mean?

Imagine we want to solve the problem A*x = b. We will do this by a permutation of the columns of A. But that does NOT involve permuting the rows of b!!!!!!!

A permutaiton of the columns of A is equivalent to a re-ordering the elements of x in that problem. That is, we expect that if x is a solution of A*x=b, and p is a column permutation of A, then we would have

A*x == A(:,p)*x(p)

There is no need to permute the elements of b, nor do you want to do so.

x2 = zeros(size(x1));
p = colamd(A);
[L,U] = lu(A(:,p)); % this builds the permutation into L and U.
x2(p) = U\(L\b);

norm(x1 - x2)
ans =
   1.1414e-13

Which works perfectly ok. The tiny difference is on the order of cond(A)*eps.

cond(full(A))*eps
ans =
   2.3173e-13

And that is why you had a problem. For some strange reason, you decided you needed to permute the elements of b. You are still wasting your time with the permutation in the first place on this tiny problem, but maybe your real problem is seriously, massively larger. Who knows?

[Balachandra Suri]

A in my case is 200000x200000 and hence it makes a huge difference doing colamd. I posted the question with 100x100 so that one can run the code and get what the problem is. Now, can you shed some insights into the problem? You can change the code to nx = 400, ny = 400 and run the code and see the difference.

I see, your updated answer addresses my issue. Thanks. That's where my math went wrong. I am also permuting b, which is incorrect. Thanks.