Delphi implementation of Aberth–Ehrlich method and precision problem

[at3s 4]

 

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.

[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.

[David Heffernan 2345]

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. 

[Anders Melander 1783]

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.

[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.

[David Heffernan 2345]

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.

[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

 

[David Heffernan 2345]

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

[Lars Fosdal 1792]

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

[Anders Melander 1783]

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.

[at3s 4]

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.

[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.

[Daniel 417]

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

[David Heffernan 2345]

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.

[Anders Melander 1783]

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 

[Anders Melander 1783]

1 hour ago, Daniel said:

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

Yes

[FPiette 383]

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).

[Lars Fosdal 1792]

Wrong floating point precision when choosing overloads
https://quality.embarcadero.com/browse/RSP-27488

Reduction of precision from integer to single rather than double
https://quality.embarcadero.com/browse/RSP-27499

Another case of silent reduction of precision
https://quality.embarcadero.com/browse/RSP-27500

[Mahdi Safsafi 225]

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. 

[David Heffernan 2345]

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.

[at3s 4]

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.

[Mahdi Safsafi 225]

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.  

[Dalija Prasnikar 1396]

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"

[at3s 4]

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?

[David Heffernan 2345]

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.