embotech / ecos

A lightweight conic solver for second-order cone programming.
GNU General Public License v3.0
468 stars 123 forks source link

Higher floating point precision than double #213

Closed DanielDoehring closed 2 years ago

DanielDoehring commented 2 years ago

I face a very ill-conditioned optimization problem where I already run for a moderate number of optimization variables into trouble. Thus, I thought about switching the floating point datatype pfloat = double too e.g. long double or, if needed, to something provided by The GNU Multiple Precision Arithmetic Library.

However, in the glbopts.h file it is explicitly stated that only double works.

Is this written as a contrast to float e.g., that (at least) double precision is required, or are any other higher precision datatypes also not possible? In the latter case, would it be possible to elaborate a bit on why that is?

adomahidi commented 2 years ago

That's an interesting direction... From an ECOS point of view there shouldn't be a fundamental issue with higher precision types, but I am not sure about the linear solver under the hood (SparseLDL). All code is there so you might give it a try and see where it fails? :)

Just curious, from which field is the optimization problem requiring more than double precision arithmetics?

DanielDoehring commented 2 years ago

Alright, I will give it a shot :)

Essentially, I am trying to optimize the coefficients $ \boldsymbol \alpha \in \mathbb R^{N+1} $ of a complex polynom $$ \sum_{i=0}^N \alphai z^i = \sum{i=0}^N \frac{\widetilde \alpha_i}{i!} z^i , \quad z \in \mathbb C $$ where $ \vert z^i \vert \sim \mathcal O (10)$. Already for $N \approx 20 $ I run into stability issues since the $z^i$ term explodes. This can be somewhat mitigated by the normalization by $i!$. However, then you have to divide the $ \tilde \alpha_i$ by $i!$ which gives you in any case $\alpha_i$ very small to zero. Of course, this problem is inherently ill-conditioned, but I want to see how far I can go with this formulation.

DanielDoehring commented 2 years ago

Switching from doubleto long double is actually not that difficult, it can be done with moderate changes. However, some adjustments have to be made, i.e., fabs becomes fabsl and printing requires usually also an "L", e.g. %f becomes %fL.

Thus, a fully type-generic version would demand some more effort.

adomahidi commented 2 years ago

Ok great to hear, so did it run with long double after those changes and did you see an improvement on numerical stability?

DanielDoehring commented 2 years ago

I called ECOS through the Julia interface where there is, unfortunately, within the Julia MOI no support for long double or other types (like Double64) but only for the standard C double.

I will now take a look at the Matlab & Python interfaces.

DanielDoehring commented 2 years ago

Based on the EiCOS project which allows fairly easy specification of constraints I forked a long double version.