Jump to content
at3s

Delphi implementation of Aberth–Ehrlich method and precision problem

Recommended Posts

 


I've implemented Aberth–Ehrlich method for finding all roots of polynomial (actually it's one of steps of eigenvalues finding).

It works very well with the degree less or equal 25.

But when degree is larger, my Delphi calculations begin to fail.

In both cases speed of calctions is acceptable (< 1 sec).

I'm quite sure it's because an accumulation of errors of very small and very large numbers during calculations.

For an algorithm needs I've developer

TComplex = record … end;

structure based on Double. The same result I got when the base data type is Extended.

I'm quite sure the problem is not an algorithm, because when I've updated TComplex to use MPArith multi precision library rather, I got exactly the same results as Octave and Matlab programs - so I assume they are correct.

But in this case time of calculation is veeery slooow (more than 20 min) because MPArith works with the memory rather.

So my questions are:

Do you know how to bypass multi precision issue in Delphi?

Do you think, if I'll port realization of calculation to another programming language (f.e. C++, C# or any other), I'll have no problem with the precision of calculations? Sure, I'll have to wrap it into dll to call within the Delphi code.

Did you already fixed a similar issue related to multi precision (maybe other more fast library)?

Any thoughts are welcome.

Share this post


Link to post
Guest
7 hours ago, at3s said:

Do you know how to bypass multi precision issue in Delphi?

Not exactly, as with Delphi i don't think it is one solution for all, but one solution per one specific problem.

 

7 hours ago, at3s said:

Do you think, if I'll port realization of calculation to another programming language (f.e. C++, C# or any other), I'll have no problem with the precision of calculations? Sure, I'll have to wrap it into dll to call within the Delphi code.

Have you tries FreePascal, and remember both FreePascal and Delphi can use the LLVM compiler which might fix your problem, but i am not sure here.

 

7 hours ago, at3s said:

Did you already fixed a similar issue related to multi precision (maybe other more fast library)?

Yes, some of them, but only for specific problem and not a general solution.

7 hours ago, at3s said:

Any thoughts are welcome.

1) Try FreePascal, then try Delphi with FMX build, FreePascal with LLVM, my recommendation then use different/better compiler like the one you validated with C++ and then link the object file, this will remove the DLL dependency.

2) As i mentioned "specific cases", so if you can paste few small test cases, the shortest the better to start with (more than one will be more accurate), where you can provide a wrong values with Delphi generated code, then someone here might help in solving it, i believe you can do it by divide and conquer, no need for full library.

Share this post


Link to post

If your goal is to find eigenvalues then you should use a dedicated library for that. There are a good few about that can handle matrices with sizes in the thousands.

 

It sounds like you are trying to find roots of the characteristic polynomial. This is not a robust approach. 

Share this post


Link to post
2 hours ago, Kas Ob. said:

Try FreePascal, then try Delphi with FMX build, FreePascal with LLVM

Why would that make a difference?

If the problem is caused by rounding errors then FP will exhibit the exact same symptoms.

Share this post


Link to post
Guest

Have a look here 

https://lemire.me/blog/2020/06/26/gcc-not-nearest/

 

Means even one compiler can generate code and produce different result with different compile settings !

So yes there is chance that FPC is doing better or worse, as OP asked for thoughts for brainstorming, so i provided some thoughts,

 

Take as example "floating-point numbers are not associative" (from that article) which not so many knows and might not need to care even, it doesn't matter that much per one value, but in the case of Aberth method it will accumulate and slide far from the right value, and i am not saying associative property is the problem, no, many things can go south.

Share this post


Link to post
4 minutes ago, Kas Ob. said:

Have a look here 

https://lemire.me/blog/2020/06/26/gcc-not-nearest/

 

Means even one compiler can generate code and produce different result with different compile settings !

So yes there is chance that FPC is doing better or worse, as OP asked for thoughts for brainstorming, so i provided some thoughts,

 

Take as example "floating-point numbers are not associative" (from that article) which not so many knows and might not need to care even, it doesn't matter that much per one value, but in the case of Aberth method it will accumulate and slide far from the right value, and i am not saying associative property is the problem, no, many things can go south.

Try other compilers isn't going to make any difference here. Other compilers will face the same problems.

 

I am reasonably confident that the correct way forward is to use the right tool for the job. And root finding of the characteristic polynomial is not the way to find eigenvalues on a finite computer.

Share this post


Link to post
Guest
13 minutes ago, David Heffernan said:

Try other compilers isn't going to make any difference here.

Most likely will not, even Wikipedia suggesting to refer to MPSolve for its better implemented method in efficiency and accuracy.

"The implementation of this method in the free software MPSolve is a reference for its efficiency and its accuracy."

 

As for Delphi and FPC, here a sample for anyone curious 

const
  xDouble   : Double = 50178230318.0;
  xExtended : Extended = 50178230318.0;
  xSingle   : Single = 50178230318.0;

  yDouble   : Double = 100000000000.0;
  yExtended : Extended = 100000000000.0;
  ySingle   : Single = 100000000000.0;

begin
  Writeln(xDouble / yDouble);
  Writeln(xExtended / yExtended);
  Writeln(xSingle / ySingle);
  Readln;
end.

Delphi Result

Quote

 5.01782303180000E-0001
 5.01782303180000E-0001
 5.01782333476502E-0001

FPC

Quote

 5.0178230318000006E-001
 5.0178230318000006E-001
 5.017823577E-01

 

Share this post


Link to post

Can it be that some of the known FP expression precision issues come into play here?

Share this post


Link to post
19 minutes ago, Kas Ob. said:

As for Delphi and FPC, here a sample for anyone curious

Yes they both suffer from floating point rounding errors. The fact that the exact error values differ doesn't make any difference. All compilers will suffer from this as the imprecision is inherent in the nature of floating point precision. Higher precision will make the error smaller but there will still be errors.

 

As David points out the only way to avoid that is to use another algorithm. If that isn't practical then one might try to implement with an arbitrary precision library.

Share this post


Link to post
1 minute ago, Anders Melander said:

Yes they both suffer from floating point rounding errors. The fact that the exact error values differ doesn't make any difference. All compilers will suffer from this as the imprecision is inherent in the nature of floating point precision. Higher precision will make the error smaller but there will still be errors.

 

As David points out the only way to avoid that is to use another algorithm. If that isn't practical then one might try to implement with an arbitrary precision library.

I'm not agree. As I mentioned above, when I used MPArith library and tested matrice 50x50, I got exactly the same result as Octave (did not test in Matlab, but there is not differences in eigenvalues results).

I also thought, I did something wrong in initial roots defining (based on Gershgorin circles).

But when I got that result with MPArith, I was convinced the problem is not in algorithm realization.

Share this post


Link to post
Guest
21 minutes ago, David Heffernan said:

@Kas Ob.It's pointless trying to argue about making the wrong solution work. 

David i am not arguing anything here.

 

OP asked and i answered, you added a though and answer, then i agreed and offered different solution, i made a point that different compilers have different result some very wrong and some close to right, and guess what some are right, from that article JavaScript is the one getting it right everytime for high precision but i didn't claim that higher precision will solve the OP problem.

 

No arguing here.

 

Just now, Anders Melander said:

Yes they both suffer from floating point rounding errors.

That wasn't the point, you claimed that they will have the same value, and i tried and regret it now to show that there is difference, now if you are going to twist it into an error is an error and no difference there, then no, i didn't claim it will solve the OP problem, and you both don't understand Aberth–Ehrlich method, it is approximation method not absolute value, this is up to OP to decide to go which way he want to go.

 

Sorry for your time, all of you.

Share this post


Link to post

Sometimes discussions suffer from the fact that participants always want to have the last word.
Several plausible causes of problems have been addressed.

Share this post


Link to post
1 hour ago, Anders Melander said:

If that isn't practical

It is practical. It is routine to solve eigen problems for systems of thousands of variables.

Share this post


Link to post
1 hour ago, Kas Ob. said:

That wasn't the point, you claimed that they will have the same value

No. I claimed that they would exhibit the same symptoms. I.e. incorrect result due to accumulation of rounding errors.

That's not twisting anything. The degree of imprecision will be the same for the same data type regardless of the compiler since we're dealing with a native data type implemented on the CPU. The exact error might differ, but that doesn't make any difference.

 

13 minutes ago, David Heffernan said:

It is routine to solve eigen problems for systems of thousands of variables.

I'll take your word for that. It's all GrΣΣk to me :classic_smile:

Share this post


Link to post
1 hour ago, Daniel said:

Sometimes discussions suffer from the fact that participants always want to have the last word.

Yes :classic_wink:

Share this post


Link to post
1 hour ago, at3s said:

I was convinced the problem is not in algorithm realization

I have not seen algorithm implementation, but one possible cause is that in your implementation, you accumulate very small floating point value and at some point the accumulated value becomes very large (or the initial value is that large). At that moment, accumulation doesn't work anymore. What very small and very large means depends on the floating point precision you used.

 

This means that how you implement the algorithm matters!

 

Compiler will not change anything, the issue comes only from the floating point precision you use (The number of significant digits).

Share this post


Link to post

My guess is that the library you're using for complexes is not well optimized as it tries to emulate complexes ! Some compilers support native complex operation such CLANG/LLVM. The wonderful Julia has a native support too. Its definitely worth checking other compilers as @Kas Ob. suggested or checking other libraries. 

Share this post


Link to post

Why does anybody think that all this guessing would be useful? Does anybody have much experience of success when guessing?

 

Solving eigen problems is a very challenging numerical problem.  Does anybody really believe that good solutions arise from guesswork?

 

I am frankly embarrassed by this thread.

  • Like 2

Share this post


Link to post

Here is part of my sources where is a main calculation:

 

    var
	  j: SmallInt;
	  nTmp, nDiv, nRootsValueOfFunction, nRootsValueOfDerivative: TComplex;
...
	function _GetValueOfFunction(ACoeff: TEquation; AValue: TComplex): TComplex; // calculate value of function F
	begin
	  Result := ACoeff.GetValueOfFunction(AValue);
	end;

	function _GetValueOfDerivative(ACoeff: TEquation; AValue: TComplex): TComplex;
	begin
	  Result := ACoeff.GetValueOfDerivative(AValue); // calculate value of derivative of function F (F')
	end;
...
	for i:= 0 to Pred(mxRoots.Cols) do
	begin      
      if not arrExcited[i] then // is not a "root" yet
      begin
        // calculate denominator
        nTmp := 0;
        for j := 0 to Pred(mxRoots.Cols) do
          if j <> i then
          begin
            nTmp := nTmp + 1 / (mxPrevRoots[0, i] - mxPrevRoots[0, j]); // mxRoots - array of root candidates; mxPrevRoots - array of "previous" root candidates
          end;
        nDiv := mxPrevRootsValueOfFunction[0, i] / mxPrevRootsValueOfDerivative[0, i]; // F/ F'
        mxRoots[0, i] := mxPrevRoots[0, i] - nDiv / (1 - nDiv * nTmp); // current calculated root value
        nRootsValueOfFunction := _GetValueOfFunction(_Self, mxRoots[0, i]); // calculate value of function F in mxRoots[0, i]
        nRootsValueOfDerivative := _GetValueOfDerivative(_Self, mxRoots[0, i]); // calculate value of derivative of function F (F') in mxRoots[0, i]
        if
          (nRootsValueOfFunction = 0)
          or
          ((mxRoots[0, i] - mxPrevRoots[0, i]).Abs <= AMistake) // detect that previous and actual roots values changed less than input error
        then
        begin
          arrExcited[i] := True; // array of Boolean; True - root is found
          Inc(iRootsCnt); // total amount of found roots
        end;
        else
        begin
          mxPrevRootsValueOfFunction[0, i] := nRootsValueOfFunction;
          mxPrevRootsValueOfDerivative[0, i] := nRootsValueOfDerivative;
        end;
	  end;
    end;

Probably someone will tell me what could be wrong here.

 

I'll try to compile my code in FPC, but I did not have experience with it and not sure if it's compilable there.

Share this post


Link to post
53 minutes ago, David Heffernan said:

Why does anybody think that all this guessing would be useful? Does anybody have much experience of success when guessing?

 

Solving eigen problems is a very challenging numerical problem.  Does anybody really believe that good solutions arise from guesswork?

 

I am frankly embarrassed by this thread.

For the same reason why anybody should think that your comment is constructive ! What does it add to OP ?

I'm going to ask you gently, please stop trolling. We just ended up a hot argument and I said I'm sorry ... Please don't let me regret.  

Share this post


Link to post

Guy is beating a dead horse with a stick.

 

Another guy comes along "Hey, your horse is dead, you need a new one"

 

"Nah... I am probably holding the stick wrong"

  • Like 2
  • Haha 2

Share this post


Link to post
10 hours ago, David Heffernan said:

If your goal is to find eigenvalues then you should use a dedicated library for that. There are a good few about that can handle matrices with sizes in the thousands.

 

It sounds like you are trying to find roots of the characteristic polynomial. This is not a robust approach.  

Can you please suggest libraries for the eigenvalues\eigenvectors finding?

Is it possible to use them in Delphi code?

Share this post


Link to post
4 hours ago, Mahdi Safsafi said:

For the same reason why anybody should think that your comment is constructive

It wasn't aimed at the OP. It was aimed at all the misleading posters who guessed and speculated. In my opinion it is those posts that are harmful to the OP's chance of progress. 

 

This topic requires knowledge of numerical methods for solving eigenproblems. Most of the posts have not lived up to that. 

Share this post


Link to post

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

×