forum.alglib.net

ALGLIB forum
It is currently Tue Jul 23, 2024 3:16 am

All times are UTC


Forum rules


1. This forum can be used for discussion of both ALGLIB-related and general numerical analysis questions
2. This forum is English-only - postings in other languages will be removed.



Post new topic Reply to topic  [ 3 posts ] 
Author Message
 Post subject: levenberg marquardt performance issue
PostPosted: Sun Feb 25, 2024 12:06 pm 
Offline

Joined: Sun Feb 25, 2024 11:37 am
Posts: 2
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
Top
 Profile  
 
 Post subject: Re: levenberg marquardt performance issue
PostPosted: Wed Feb 28, 2024 12:32 pm 
Offline
Site Admin

Joined: Fri May 07, 2010 7:06 am
Posts: 921
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.


Top
 Profile  
 
 Post subject: Re: levenberg marquardt performance issue
PostPosted: Wed Feb 28, 2024 9:04 pm 
Offline

Joined: Sun Feb 25, 2024 11:37 am
Posts: 2
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


Top
 Profile  
 
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 3 posts ] 

All times are UTC


Who is online

Users browsing this forum: No registered users and 2 guests


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot post attachments in this forum

Search for:
Jump to:  
cron
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group