forum.alglib.net
http://forum.alglib.net/

levenberg marquardt performance issue
http://forum.alglib.net/viewtopic.php?f=2&t=4593
Page 1 of 1

Author:  lesauxvi [ Sun Feb 25, 2024 12:06 pm ]
Post subject:  levenberg marquardt performance issue

Dear all,

I am searching for alternatives to mpfit for a C++ program. The optimization subpackage of alglib seems to be a good candidate. Before going any further, I wrote a small example where the same problem is solved using both libraries. I managed to get the same result, but alglib is way slower (factor 8). Before jumping to conclusions, I want to make sure that I haven't missed any important things that could impact on performance. The code is compiled with cmake, and I did nothing specific in terms of compilation flag, except a release mode -O3 flag.

Could an alglib expert see what I'm doing wrong? I have also added a full example in attachement. I had to remove the alglib sources to keep the file under 1Mo, but you just have to copy/paste the sources in the alglib folder and it should work.

If the commercial edition brings a performance boost, I would definitely consider bying it.

Here is the code:

Code:
#include <algorithm>
#include <cmath>
#include <chrono>
#include <iostream>
#include <vector>
// alglib
#include "src/ap.h"
#include "src/stdafx.h"
#include "src/optimization.h"

struct vars_struct
{
    double *x;
    double *y;
};

void planck_alglib(const alglib::real_1d_array &p, alglib::real_1d_array &vy, void *ptr)
{
    struct vars_struct *v = (struct vars_struct *) ptr;
    for (int i=0; i<vy.length(); i++)
        vy[i] = ((p[0] / (std::exp(p[1] / v->x[i]) - 1) + p[2]) - v->y[i]);
}


void  planck_alglib_jac(const alglib::real_1d_array &p, alglib::real_1d_array &vy, alglib::real_2d_array &jac, void *ptr)
{
    struct vars_struct *v = (struct vars_struct *) ptr;
    double B_X;
    for (int i=0; i<vy.length(); i++) {
        B_X = std::exp(p[1] / v->x[i]);
        vy[i] = (p[0] / (B_X - 1) + p[2]) - v->y[i];
        jac[i][0] = 1 / (B_X - 1);
        jac[i][1] = - p[0] * B_X / ((B_X - 1) * (B_X - 1) * v->x[i]);
        jac[i][2] = 1;
    }
}



int main() {
        // 50 data points
        std::vector<double> T = {
        -5.38739, -3.7718, -2.23958, -0.78166, 0.60961, 1.94071, 3.21719, 4.44385,
        5.62487, 6.76392, 7.86421, 8.9286, 9.95963, 10.95957, 11.93047, 12.87417,
        13.79233, 14.68647, 15.55796, 16.40807, 17.23796, 18.04867, 18.84118, 19.61641,
        20.37518, 21.11827, 21.84639, 22.56021, 23.26036, 23.94741, 24.62192, 25.2844,
        25.93531, 26.57511, 27.20423, 27.82306, 28.43198, 29.03135, 29.62149, 30.20273,
        30.77538, 31.3397, 31.89598, 32.44447, 32.98542, 33.51906, 34.04562, 34.5653,
        35.0783, 35.58483};
    std::transform(T.begin(), T.end(), T.begin(), [&](double& value){ return value + 273.15; });
    std::vector<double> DV = {
        3762, 3910, 4058, 4206, 4354, 4502, 4650, 4798, 4946,5094, 5242, 5390, 5538,
        5686, 5834, 5982, 6130, 6278, 6426, 6574, 6722, 6870, 7018, 7166, 7314, 7462,
        7610, 7758, 7906, 8054, 8202, 8350, 8498, 8646, 8794, 8942, 9090, 9238, 9386,
        9534, 9682, 9830, 9978, 10126, 10274, 10422, 10570, 10718, 10866, 11014};

    struct vars_struct v;
    v.x = T.data();
    v.y = DV.data();
   
    alglib::real_1d_array s = "[100000, 1, 1]";
    alglib::real_1d_array p0 = "[100000000, 1000, 1000]";
    alglib::minlmstate state;
    double epsx = 0.000001;
    alglib::ae_int_t maxits = 200;
    alglib::minlmreport rep;
    alglib::minlmcreatevj(3, 50, p0, state);
    alglib::minlmsetcond(state, epsx, maxits);
    alglib::minlmsetscale(state, s);
    begin = std::chrono::steady_clock::now();
    // we simulate here various optimization, to check if minlmrestartfrom can help us gain some time
    for (int n = 0; n < 100; n++) {
        p0[0] = 100000000;
        p0[1] = 1000;
        p0[2] = 1000;
        alglib::minlmrestartfrom(state, p0);
        alglib::minlmoptimize(state, planck_alglib, planck_alglib_jac, nullptr, (void*)&v);
        minlmresults(state, p0, rep);
    }
    end = std::chrono::steady_clock::now();

    std::cout << "R: " << p0[0] << std::endl;   // R: 1.476e+08
    std::cout << "B: " << p0[1] << std::endl;   // B: 2983.1
    std::cout << "O: " << p0[2] << std::endl;  // O: 1620.77
    std::cout << "duration : " << std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count() << " ms" << std::endl;

    return 0;
}


I woud really appreciate if someone could help me understand what is going on. Many thanks in advance

Vincent

Attachments:
ful_example.zip [61.19 KiB]
Downloaded 163 times

Author:  Sergey.Bochkanov [ Wed Feb 28, 2024 12:32 pm ]
Post subject:  Re: levenberg marquardt performance issue

Hi! There is nothing wrong with your code, and ALGLIB is ok too. Two reasons for the difference in performance are:

* ALGLIB spends a lot of time in the internal integrity checks, and on such small problems it constitutes a non-negligible fraction of the running time. MPFIT is quite fast on small and easy problems like yours, but it has no integrity checks like code for handling the model breaking down (returning infinities) or other problematic cases.

* MPFIT needs 4x less iterations than ALGLIB because it uses diag(J'J) as a Levenberg-Marquardt scaling factor, and ALGLIB uses identity as a scaling factor. The former (MPFIT) gives much faster convergence, but the latter (ALGLIB) is more robust. diag(J'J) as a step scaling sometimes fails on difficult problems, and basically every developer who implements LM solver has to choose between performance and robustness here.

Author:  lesauxvi [ Wed Feb 28, 2024 9:04 pm ]
Post subject:  Re: levenberg marquardt performance issue

Dear Sergey,

Many thanks for your valuable feedback and for taking the time to reply. I totally agree on the fact that my problem is very simple (no to say trivial), but I have to repeat it over a lot of time (1 million times usually) with the same type of data (same number of points). I initially thought that the restart functionnality of alglib could help me gain some time (this functionnality is not in mpfit).

I will keep mpfit for now but I keep in mind alglib for other purposes.

Vincent

Page 1 of 1 All times are UTC
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group
http://www.phpbb.com/