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 |
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 |
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 |
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 |
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 |
Free forum by Nabble | Edit this page |