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,
> 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
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.
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
On May 15, 12:07 am, Alois Steindl <Alois.Stei...@tuwien.ac.at> wrote:
> 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
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,
<goodchild.tre...@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).
> 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,
<goodchild.tre...@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.
> > 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
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
>> 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
Spellucci <spellu...@fb04373.mathematik.tu-darmstadt.de> wrote: > In article <200520082258269374%cbarron...@adelphia.net>, > Carl Barron <cbarron...@adelphia.net> writes: > >In article > ><416e63ba-a9a7-44e7-917c-5882e78c9...@34g2000hsh.googlegroups.com>, > ><goodchild.tre...@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
> 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
