Solving a large, sparse, nearly symmetric system

[goodchil...@gmail.com]

Hi,

I'm solving a very large, real, sparse system that has the following
structure:

[ A B^T ]
[ B 0 ]

A is banded, nonzero along the main diagonal, and the off-diagonals
are antisymmetric. In other words, ignoring the main diagonal, A = -
A^T. Additionally, A is diagonally dominant. Finally, B B^T is
symmetric positive definite.

I would like to avoid methods that potentially require very dense
factors since I'm solving fairly large systems. Ideally I'd like to
use a method for which there's an existing, efficient, freely
available implementation, or one that's reasonably easy to implement
efficiently. The only thing I know of that satisfies all my
constraints is biconjugate gradient or some simple iterative solver
like Jacobi or Gauss-Seidel. As you can probably tell, I'm not an
expert in this field, so all suggestions are greatly appreciated!
Thanks,

Trevor

[goodchil...@gmail.com]

One other potentially useful piece of information: I'm solving a long
sequence of systems of this form wherein A changes but B is always the
same.

-Trevor

[Alois Steindl]

Hello,
you stated a quite similar question a short time ago.
Usually my intention to help decreases significantly, if someone
simply ignores the hints, which I gave to the first question and just
restates the same question again.
What is wrong about using BEMW?
Alois

--
Alois Steindl, Tel.: +43 (1) 58801 / 32558
Inst. for Mechanics and Mechatronics Fax.: +43 (1) 58801 / 32598
Vienna University of Technology, A-1040 Wiedner Hauptstr. 8-10

[goodchil...@gmail.com]

Hi Alois,

Sorry if I seem to have ignored your suggestions -- actually, the
question I asked earlier was about an entirely different problem (one
involving a dense column/dense row, and without any of the special
structure of this problem). I have also tried to apply your
suggestion from the earlier problem of writing out the equations
resulting from the block form and solving a series of two systems.
This approach led me to look at null space methods which seemed to
offer some performance advantage, but didn't really seem to be
exploiting much of the special structure I have in this latest problem
(and added quite a bit of complexity in determining the null space
matrix). My hope was that since the problem seems to be very simple
and have a lot of special structure I would be able to use some
``tricks'' that are more specific to this system than simply running
BEMW. Again, this is not my area of expertise, so I'm sorry if this
is a strange request, or if I'm missing an obvious connection between
this problem and my previous one.
Thanks for your help,

Trevor

[Carl Barron]

<goodchil...@gmail.com> wrote:

> On May 14, 6:14 pm, goodchild.tre...@gmail.com wrote:
> > Hi,
> >
> > I'm solving a very large, real, sparse system that has the following
> > structure:
> >
> > [ A B^T ]
> > [ B 0 ]
> >
> > A is banded, nonzero along the main diagonal, and the off-diagonals
> > are antisymmetric. In other words, ignoring the main diagonal, A = -
>

> > Trevor
>
> One other potentially useful piece of information: I'm solving a long
> sequence of systems of this form wherein A changes but B is always the
> same.
>
> -Trevor

Assuming A is m x m and B is n x m m much biger than n. One can
solve two sets of n equations and n unknowns both with with the same
left hand size BB^T which is n x n and psd.

if we partition the solution vector into (x,y)^T so the following makes
sense you have:
Ax + By = b and Bx = d.

let B^Tz = x then BB^Tz = d since BB^T is psd a choleski
factorization exists. LL^T = BB^T This is done once.

so now ew find z and x = B^Tz

Now we have B^Ty = b - Ax.

BB^Ty = B(b-Ax)
which finds y by uising LL^T above.

done once
find BB^T O(n^2m)
factor BB^T O(n^3)
done once for each changind d
forward / back solve LL^Tz = d O(n^2)
find x = B^Tz O(nm)
for each changing A and/or b.
find b - Ax O(mp) p = bandwidth.
find B(b-Ax) O(mn)
find y O(n^2) forward /backward solutions.

so if n is smallish say ~ 100 this involves a dense solver of a 100 x
109 system. not too taxing as the oprerations invoving m sized variables
are just vector multiplies.

operations of O(n^3) are done once. no operations O(m^3).

[goodchil...@gmail.com]

On May 16, 9:17 pm, cbarr...@ix.netcom.com (Carl Barron) wrote:

> Assuming A is m x m and B is n x m m much biger than n. One can
> solve two sets of n equations and n unknowns both with with the same
> left hand size BB^T which is n x n and psd.
>
> if we partition the solution vector into (x,y)^T so the following makes
> sense you have:
> Ax + By = b and Bx = d.
>
> let B^Tz = x then BB^Tz = d since BB^T is psd a choleski
> factorization exists. LL^T = BB^T This is done once.
>
> so now ew find z and x = B^Tz
>
> Now we have B^Ty = b - Ax.
>
> BB^Ty = B(b-Ax)
> which finds y by uising LL^T above.

Hi Carl,

Thanks for your solution. Unfortunately, you make one assumption
about the problem that I cannot make - namely that x is in the range
of B^T (i.e., "let B^Tz = x"). I happen to know in my problem that
this assumption is not valid, and it would be quite strange if x did
not depend at all on A in my particular problem.
Thanks again for the help,

Trevor

[Carl Barron]

In article
<416e63ba-a9a7-44e7...@34g2000hsh.googlegroups.com>,
<goodchil...@gmail.com> wrote:

The only thing I can come up with requres n+1 m dimensional band
solvers and one n x n dense system but does not take advantage of B
being fixed.
| I 0 | |A B^T | = |A B^T| CA = B
|-C I | |B 0 | | 0 D | D = -CB^T

|A B^T | |x1| = | b|
|0 D | |x2| | d - Cb|

x = |x1|
|x2|

but I really dont like it.

[Peter Spellucci]

In article <200520082258269374%cbarr...@adelphia.net>,

Carl Barron <cbarr...@adelphia.net> writes:
>In article
><416e63ba-a9a7-44e7...@34g2000hsh.googlegroups.com>,
><goodchil...@gmail.com> wrote:
>
>> On May 16, 9:17 pm, cbarr...@ix.netcom.com (Carl Barron) wrote:
>>
>> > Assuming A is m x m and B is n x m m much biger than n. One can
>> > solve two sets of n equations and n unknowns both with with the same
>> > left hand size BB^T which is n x n and psd.
>> >

snip

> The only thing I can come up with requres n+1 m dimensional band
>solvers and one n x n dense system but does not take advantage of B
>being fixed.
>| I 0 | |A B^T | = |A B^T| CA = B
>|-C I | |B 0 | | 0 D | D = -CB^T
>
>|A B^T | |x1| = | b|
>|0 D | |x2| | d - Cb|
>
>x = |x1|
> |x2|
>
>but I really dont like it.

since you have no symmetry in A (but strict diagonal dominance if I remember
right, you even cannot apply the special software for KKT systems (there is
now a vast literature about this)) but you could do an ordinary sparse
LU decomposition (without pivoting first, evenyually after rescaling B ,
which is of no promlem in your case). then you could use this one as preconditioner
for an iterative solver if A has changed (a little)
e.g. bicg-stab. the usual way to exploit the properties you mention would be:
use a sparse QR decomposition of B, transform the system using this
(this fixes x partly) and then to resolve "for the rest". For A you could
use a A= LU decomposition without pivoting. but then
it remains so to solve a subsystem of the form

[O,I_{m-n}]QLUQ'[O;I_{n-m}] tilde x_2 = tilde b_2

and although you could have Q represented in sparse Givens or Householder
factor and L and U banded and sparse the matrix on the left would usually
be dense. hence this approach would make sense only in an iterative solver
scheme but this will work only if the spectrum of A is well behaved and
which does not need to compute the matrix explicitly.
you have the relation
Q x = tilde x
and tilde x_2 are the last n-m components of it. maybe you can invent
a special iterative scheme for this ?
hth
peter

[Carl Barron]

In article <g110ks$c1j$1...@fb04373.mathematik.tu-darmstadt.de>, Peter
Spellucci <spel...@fb04373.mathematik.tu-darmstadt.de> wrote:

The op might be interested in this link as it seems very related to
the problem:
www.precond07.enseeiht.fr/Talks/sameh/sameh.pdf
google sameh fluid dynamics will find it also.