about sovling equation Ax=b

[sl]

I need to solve Ax=b, while A is a large sparse matrix. some
information about A

A : (n+k)×n , k usually is 3, or 4 etc.
rank(A)=n
and I need to sovle Ax=b many times for different RHS b

anyone have suggestion for me? which tool is suitable for me.

I think the best way is first factor A the solve the equation for
different b/

how to do it?

[Mark Hoemmen]

1. Since A has more rows than columns, you may not be able to solve
Ax=b exactly. Are you interested in a least-squares solution?

2. Since A is large and sparse but does not have full column rank,
you'll need a special kind of factorization (ordinary sparse QR won't
suffice). I suspect (but am not sure, as this is not my area of
expertise) that you will need a sparse rank-revealing factorization.

3. Can you tell us more about your application? This might help us
solve your problem more efficiently.

mfh

[sl]

Thank you , Mark
I do research in shape deformation. And need to solve Ax=b for A
large sparse matrix.

[Evgenii Rudnyi]

sl schrieb:

> Thank you , Mark
> I do research in shape deformation. And need to solve Ax=b for A
> large sparse matrix.

You have to define your mathematical problem. Usually the equation Ax =
B defines a system of linear equations where the number of equations is
equal to the number of unknowns. In this case just use one of the sparse
solvers, you may find some links here

http://MatrixProgramming.com

If however the number of equations is not equal to the number of
unknowns, then one should define what is the goal.

[Mark Hoemmen]

On Sat, Aug 29, 2009 at 08:39, Evgenii Rudnyi<use...@rudnyi.ru> wrote:
> If however the number of equations is not equal to the number of
> unknowns, then one should define what is the goal.

It's important also to point out that since your matrix is
rank-deficient, it depends on your problem whether you expect to be
able to solve an equation Ax=b "exactly," or whether you really need
to pick x to minimize some measure(s) or error. The latter is usually
what you want.

mfh

[sl]

<<or whether you really need to pick x to minimize some measure(s) or
error. The latter is usually
what you want.

You are right . because the number of equation is bigger than the
number of unkonwn
the solution is at the sense of least square.

[Mark Hoemmen]

On Sat, Aug 29, 2009 at 19:06, sl<hit...@qq.com> wrote:
> You are right . because the number of equation is bigger than the

> number of unknowns, the solution is at the sense of least squares.

Since the matrix A does not have full column rank, generally one also
imposes the condition that the solution x should have minimum 2-norm.
Is that appropriate for your application?

mfh

[Perry]

I would like to activate this question (not fully resolved here) as I need to solve a very similar problem.

The problem can be restated as:

"This is equivalent to solving a sparse linear system Ax = b in least squares sense. Thus x can be solved from the normal equations A^TA x = A^Tb."

There is one step missing here: from sparse matrix A (size of (n+k)*n, n*n part being symmetric), how to compute A^TA and A^Tb and save the results as TAUCS compatible and ultimately solve the linear system using TAUCS.


On Sunday, 30 August 2009 13:21:42 UTC-4, Mark H. wrote:

[Evgenii Rudnyi]

An example of sparse matrix multiply is here

http://matrixprogramming.com/files/code/Sparse/

Yet, it may not be efficient.

Evgenii

On 18.07.2012 00:52 Perry said the following:

[Perry]

Thanks Evgenii.

Your "class SparseMatrix" assumes the matrix to be a square one. Therefore I have the problem of computing the multiplications of A^TA and A^Tb, where A is an (n+k)*n sparse matrix and b is an (n+k) dense vector. Please advise.

[Evgenii Rudnyi]

On 20.07.2012 01:09 Perry said the following:

> Thanks Evgenii.
>
> Your "class SparseMatrix" assumes the matrix to be a square one.
> Therefore I have the problem of computing the multiplications of A^TA
> and A^Tb, where A is an (n+k)*n sparse matrix and b is an (n+k) dense
> vector. Please advise.
>

Correct. It is necessary to modify it. The code was just an example, on
how one could do it, I do not have a ready solution for you, sorry.

Evgenii

[Perry]

Hi Evgenii,

I am trying to modify the code you have provided. I thus have questions on the usages of "swap()" function in part of your code:

"
void swap(SparseMatrix &y)
    {
        vector< map<int, double> >::swap(y);      // why swap only takes in one input?
        std::swap(m, y.m);                                    // don't understand the purpose of this line either
    }
    void clearMemory()
    {
        SparseMatrix empty;
        swap(empty);                                           // again why only one input?
    }
"
The above was defined in the class of SparseMatrix; however, it is not used in the example code of SparseMatrix. Could you explain what it is defined for? According to the grammar of swap function, it usually takes two inputs as parameters but I see here you only used one. Please explain.

By reading into your code, I am getting to understand the way how sparse matrix is organized and used for computation. Therefore this makes it possible for me to modify your code for my problem. The above questions may be trivial (or not) but I would like to get it clear in case I missed anything in your code.

Keeping the sparse matrix by rows in your code is great as it facilitates the computation of C = AA'. If I were to compute C=A'A, it may be better for me to first transpose A and then call the function of C=AA'. This should be easier than keeping the sparse matrix by cols and then defining C=A'A. Please verify.

Again thanks for leading me through the sparse matrix computation.

Perry

[Evgenii Rudnyi]

On 24.07.2012 18:03 Perry said the following:

> Hi Evgenii,
>
> I am trying to modify the code you have provided. I thus have
> questions on the usages of "swap()" function in part of your code:

The function swap makes the change between two objects indeed. Yet, the
function belongs to an object and this is the first argument (this).

For example

SparseMatrix A;
SparseMatrix B;

A.swap(B);

The last operation changes the content between A and B. It is the same as

B.swap(A);

> " void swap(SparseMatrix&y) { vector< map<int, double> >::swap(y);

> // why swap only takes in one input? std::swap(m, y.m);

The above is different from the case when one program a function that do
not belong to the class. This could be also possible to develop.

> // don't understand the purpose of this line either }


>void clearMemory()
{
SparseMatrix empty;
swap(empty);
}

Here it is same. You have an object

SparseMatrix A;

A.clearMemory();

will remove all arrays associated with A.

Evgenii