QuantLib › quantlib-dev

BiCGstab

Classic

List

Threaded

5 messages Options

Peter Caspers-4

BiCGstab

Hi Ralph, Klaus,

I am using your biconjugate gradient stabilized solver and experience
some issues in certain cases. A small test case would be the system
Ax=b with

A=
| 1 0 0 0 0 0 0 0 0 |
| 0 1 0 0 0 0 0 0 0 |
| 0 0 1 0 0 0 0 0 0 |
| 0 0 0 1 0 0 0 0 0 |
| 0 0 0 0 1 0 0 0 0 |
| 0 0 -0.5 0 0 1 0 0 -0.5 |
| 0 0 0 -0.5 0 0 1 -0.5 0 |
| 0 0 0 0 0 0 0 1 0 |
| 0 0 0 0 0 -0.5 0 -0.5 1 |

b = [ 5; 5; 5; 5; 5; 0; 0; 5; 0 ]

With the guess [5,5,...,5], which is the exact solution, the solver
throws an exception "could not converge". This is because "error" is
initialized to Null<Real>() in L52, and the iteration loop is exited
at L58 already (because rho=0.0), so the Null<Real>() value is checked
in L90. I think this can be fixed by simply changing L52 to

Real error=norm2(r)/bnorm2;

With the guess [0,0,....,0] on the other hand, the same exception is
thrown. We have

x = [ 5; 5; 5; 5; 5; 2.77778; 5.55556; 5; 2.77778 ]
r = [ 0; 0; 0; 0; 0; 1.11111; -0.555556; 0; 1.11111 ]
rTld = [ 5; 5; 5; 5; 5; 0; 0; 5; 0 ]

at L86 after the first iteration, so rho = 0.0 again, but this time
exiting (at L57) with a huge error of course. Is that expected ? While
I have not looked at the code very closely, a quick observation is
that rTld is never updated during the iterations, which doesn't feel
right intuitively.

Can you help?

Thanks
Peter

------------------------------------------------------------------------------
_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

Ralph Schreyer

Re: BiCGstab

Hi Peter,

I've just had a quick look ... The proper initialization of "error" seems to be missing and changing L52 to

Real error=norm2(r)/bnorm2;

as you suggested would be the correct bugfix.

On the other issue, the intermediate results rho == 0.0 or omega == 0.0 are so-called breakdown cases for the algorithm. To cope with that, one would have to combine the BiCG with other methods ... In the BiCGstab algorithm, rTld is not updated indeed.

Best regards

Ralph

Gesendet: Mittwoch, 18. November 2015 um 19:53 Uhr
Von: "Peter Caspers" <[hidden email]>
An: "QuantLib Mailing Lists" <[hidden email]>
Betreff: [Quantlib-dev] BiCGstab

------------------------------------------------------------------------------

_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

Peter Caspers-4

Re: BiCGstab

Hi Ralph,

thanks a lot. I will send Luigi a PR with the fix for the error initialization.

Best regards
Peter

On 19 November 2015 at 15:03, <[hidden email]> wrote:

> Hi Peter,
>
> I've just had a quick look ... The proper initialization of "error" seems to
> be missing and changing L52 to
>
> Real error=norm2(r)/bnorm2;
>
> as you suggested would be the correct bugfix.
>
> On the other issue, the intermediate results rho == 0.0 or omega == 0.0 are
> so-called breakdown cases for the algorithm. To cope with that, one would
> have to combine the BiCG with other methods ... In the BiCGstab algorithm,
> rTld is not updated indeed.
>
>
> Best regards
>
> Ralph
>
>
>
>
> Gesendet: Mittwoch, 18. November 2015 um 19:53 Uhr
> Von: "Peter Caspers" <[hidden email]>
> An: "QuantLib Mailing Lists" <[hidden email]>
> Betreff: [Quantlib-dev] BiCGstab
> Hi Ralph, Klaus,
>
> I am using your biconjugate gradient stabilized solver and experience
> some issues in certain cases. A small test case would be the system
> Ax=b with
>
> A=
> | 1 0 0 0 0 0 0 0 0 |
> | 0 1 0 0 0 0 0 0 0 |
> | 0 0 1 0 0 0 0 0 0 |
> | 0 0 0 1 0 0 0 0 0 |
> | 0 0 0 0 1 0 0 0 0 |
> | 0 0 -0.5 0 0 1 0 0 -0.5 |
> | 0 0 0 -0.5 0 0 1 -0.5 0 |
> | 0 0 0 0 0 0 0 1 0 |
> | 0 0 0 0 0 -0.5 0 -0.5 1 |
>
> b = [ 5; 5; 5; 5; 5; 0; 0; 5; 0 ]
>
> With the guess [5,5,...,5], which is the exact solution, the solver
> throws an exception "could not converge". This is because "error" is
> initialized to Null<Real>() in L52, and the iteration loop is exited
> at L58 already (because rho=0.0), so the Null<Real>() value is checked
> in L90. I think this can be fixed by simply changing L52 to
>
> Real error=norm2(r)/bnorm2;
>
> With the guess [0,0,....,0] on the other hand, the same exception is
> thrown. We have
>
> x = [ 5; 5; 5; 5; 5; 2.77778; 5.55556; 5; 2.77778 ]
> r = [ 0; 0; 0; 0; 0; 1.11111; -0.555556; 0; 1.11111 ]
> rTld = [ 5; 5; 5; 5; 5; 0; 0; 5; 0 ]
>
> at L86 after the first iteration, so rho = 0.0 again, but this time
> exiting (at L57) with a huge error of course. Is that expected ? While
> I have not looked at the code very closely, a quick observation is
> that rTld is never updated during the iterations, which doesn't feel
> right intuitively.
>
> Can you help?
>
> Thanks
> Peter
>
> ------------------------------------------------------------------------------
> _______________________________________________
> QuantLib-dev mailing list
> [hidden email]
> https://lists.sourceforge.net/lists/listinfo/quantlib-dev

------------------------------------------------------------------------------
_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

Klaus Spanderen-2

Re: BiCGstab

Hi Peter,

rTld does not get updated, which makes the choice of the guess even more
important. If b_i and (A*guess)_i are zero then BICGSTAB is likely to fail if
(A^(-1)*b)_i is not zero.

MatLab also fails to solve your second example using MatLab's own BICGSTAB
implementation.

regards
Klaus

On Thursday, November 19, 2015 03:55:54 PM Peter Caspers wrote:

> Hi Ralph,
>
> thanks a lot. I will send Luigi a PR with the fix for the error
> initialization.
>
> Best regards
> Peter
>
> On 19 November 2015 at 15:03, <[hidden email]> wrote:
> > Hi Peter,
> >
> > I've just had a quick look ... The proper initialization of "error" seems
> > to be missing and changing L52 to
> >
> > Real error=norm2(r)/bnorm2;
> >
> > as you suggested would be the correct bugfix.
> >
> > On the other issue, the intermediate results rho == 0.0 or omega == 0.0
> > are
> > so-called breakdown cases for the algorithm. To cope with that, one would
> > have to combine the BiCG with other methods ... In the BiCGstab algorithm,
> > rTld is not updated indeed.
> >
> >
> > Best regards
> >
> > Ralph
> >
> >
> >
> >
> > Gesendet: Mittwoch, 18. November 2015 um 19:53 Uhr
> > Von: "Peter Caspers" <[hidden email]>
> > An: "QuantLib Mailing Lists" <[hidden email]>
> > Betreff: [Quantlib-dev] BiCGstab
> > Hi Ralph, Klaus,
> >
> > I am using your biconjugate gradient stabilized solver and experience
> > some issues in certain cases. A small test case would be the system
> > Ax=b with
> >
> > A=
> >
> > | 1 0 0 0 0 0 0 0 0 |
> > | 0 1 0 0 0 0 0 0 0 |
> > | 0 0 1 0 0 0 0 0 0 |
> > | 0 0 0 1 0 0 0 0 0 |
> > | 0 0 0 0 1 0 0 0 0 |
> > | 0 0 -0.5 0 0 1 0 0 -0.5 |
> > | 0 0 0 -0.5 0 0 1 -0.5 0 |
> > | 0 0 0 0 0 0 0 1 0 |
> > | 0 0 0 0 0 -0.5 0 -0.5 1 |
> >
> > b = [ 5; 5; 5; 5; 5; 0; 0; 5; 0 ]
> >
> > With the guess [5,5,...,5], which is the exact solution, the solver
> > throws an exception "could not converge". This is because "error" is
> > initialized to Null<Real>() in L52, and the iteration loop is exited
> > at L58 already (because rho=0.0), so the Null<Real>() value is checked
> > in L90. I think this can be fixed by simply changing L52 to
> >
> > Real error=norm2(r)/bnorm2;
> >
> > With the guess [0,0,....,0] on the other hand, the same exception is
> > thrown. We have
> >
> > x = [ 5; 5; 5; 5; 5; 2.77778; 5.55556; 5; 2.77778 ]
> > r = [ 0; 0; 0; 0; 0; 1.11111; -0.555556; 0; 1.11111 ]
> > rTld = [ 5; 5; 5; 5; 5; 0; 0; 5; 0 ]
> >
> > at L86 after the first iteration, so rho = 0.0 again, but this time
> > exiting (at L57) with a huge error of course. Is that expected ? While
> > I have not looked at the code very closely, a quick observation is
> > that rTld is never updated during the iterations, which doesn't feel
> > right intuitively.
> >
> > Can you help?
> >
> > Thanks
> > Peter
> >
> > --------------------------------------------------------------------------
> > ---- _______________________________________________
> > QuantLib-dev mailing list
> > [hidden email]
> > https://lists.sourceforge.net/lists/listinfo/quantlib-dev
>
> ----------------------------------------------------------------------------
> -- _______________________________________________
> QuantLib-dev mailing list
> [hidden email]
> https://lists.sourceforge.net/lists/listinfo/quantlib-dev

------------------------------------------------------------------------------
_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

Peter Caspers-4

Re: BiCGstab

Hi Klaus,

thank you. Actually the breakdown test case was born because I forgot
to provide the guess in the solve(...) method (and then 0 is taken).
With the guess provided (and the fix discussed earlier) I did not run
into any problems up to now, also in bigger and more realistic / less
trivial cases for A and b.

In the Numerical Recipes (2.7) I found this statement

"While there is no guarantee that this whole procedure will not break
down or become
unstable for general A, in practice this is rare."

So everything's fine, sorry for my slight disbelief in the
implementation, it was based on ignorance :-)

Thanks again for your help
Peter

On 19 November 2015 at 21:31, Klaus Spanderen <[hidden email]> wrote:

> Hi Peter,
>
> rTld does not get updated, which makes the choice of the guess even more
> important. If b_i and (A*guess)_i are zero then BICGSTAB is likely to fail if
> (A^(-1)*b)_i is not zero.
>
> MatLab also fails to solve your second example using MatLab's own BICGSTAB
> implementation.
>
> regards
> Klaus
>
> On Thursday, November 19, 2015 03:55:54 PM Peter Caspers wrote:
>> Hi Ralph,
>>
>> thanks a lot. I will send Luigi a PR with the fix for the error
>> initialization.
>>
>> Best regards
>> Peter
>>
>> On 19 November 2015 at 15:03, <[hidden email]> wrote:
>> > Hi Peter,
>> >
>> > I've just had a quick look ... The proper initialization of "error" seems
>> > to be missing and changing L52 to
>> >
>> > Real error=norm2(r)/bnorm2;
>> >
>> > as you suggested would be the correct bugfix.
>> >
>> > On the other issue, the intermediate results rho == 0.0 or omega == 0.0
>> > are
>> > so-called breakdown cases for the algorithm. To cope with that, one would
>> > have to combine the BiCG with other methods ... In the BiCGstab algorithm,
>> > rTld is not updated indeed.
>> >
>> >
>> > Best regards
>> >
>> > Ralph
>> >
>> >
>> >
>> >
>> > Gesendet: Mittwoch, 18. November 2015 um 19:53 Uhr
>> > Von: "Peter Caspers" <[hidden email]>
>> > An: "QuantLib Mailing Lists" <[hidden email]>
>> > Betreff: [Quantlib-dev] BiCGstab
>> > Hi Ralph, Klaus,
>> >
>> > I am using your biconjugate gradient stabilized solver and experience
>> > some issues in certain cases. A small test case would be the system
>> > Ax=b with
>> >
>> > A=
>> >
>> > | 1 0 0 0 0 0 0 0 0 |
>> > | 0 1 0 0 0 0 0 0 0 |
>> > | 0 0 1 0 0 0 0 0 0 |
>> > | 0 0 0 1 0 0 0 0 0 |
>> > | 0 0 0 0 1 0 0 0 0 |
>> > | 0 0 -0.5 0 0 1 0 0 -0.5 |
>> > | 0 0 0 -0.5 0 0 1 -0.5 0 |
>> > | 0 0 0 0 0 0 0 1 0 |
>> > | 0 0 0 0 0 -0.5 0 -0.5 1 |
>> >
>> > b = [ 5; 5; 5; 5; 5; 0; 0; 5; 0 ]
>> >
>> > With the guess [5,5,...,5], which is the exact solution, the solver
>> > throws an exception "could not converge". This is because "error" is
>> > initialized to Null<Real>() in L52, and the iteration loop is exited
>> > at L58 already (because rho=0.0), so the Null<Real>() value is checked
>> > in L90. I think this can be fixed by simply changing L52 to
>> >
>> > Real error=norm2(r)/bnorm2;
>> >
>> > With the guess [0,0,....,0] on the other hand, the same exception is
>> > thrown. We have
>> >
>> > x = [ 5; 5; 5; 5; 5; 2.77778; 5.55556; 5; 2.77778 ]
>> > r = [ 0; 0; 0; 0; 0; 1.11111; -0.555556; 0; 1.11111 ]
>> > rTld = [ 5; 5; 5; 5; 5; 0; 0; 5; 0 ]
>> >
>> > at L86 after the first iteration, so rho = 0.0 again, but this time
>> > exiting (at L57) with a huge error of course. Is that expected ? While
>> > I have not looked at the code very closely, a quick observation is
>> > that rTld is never updated during the iterations, which doesn't feel
>> > right intuitively.
>> >
>> > Can you help?
>> >
>> > Thanks
>> > Peter
>> >
>> > --------------------------------------------------------------------------
>> > ---- _______________________________________________
>> > QuantLib-dev mailing list
>> > [hidden email]
>> > https://lists.sourceforge.net/lists/listinfo/quantlib-dev
>>
>> ----------------------------------------------------------------------------
>> -- _______________________________________________
>> QuantLib-dev mailing list
>> [hidden email]
>> https://lists.sourceforge.net/lists/listinfo/quantlib-dev
>

------------------------------------------------------------------------------
_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev