trying to set up Gauss-Newton method for fitting a power law curve

[American Citizen <website.reader3@gmail.com>]

Hello:

The count of triads of squares = n (a given number) seems to follow a 
power law of n = 5.9 * count ^ 1.92 more or less.

Here is my attempt at implementing that algorithm, but it is running 
very slowly for 80,942 elements [ count, n ] vector.

Is this about right?

Randall

/* Function to calculate the residuals based on current parameters a and 
b */
residuals(a, b, data) = {
    my(res, n, count, i);
    res = 0; \\ Initialize the residual sum
    for (i = 1, length(data),
        n = data[i][2]; \\ Get n
        count = data[i][1]; \\ Get count
        res += (count - a * n^b)^2; \\ Calculate the squared residual
    );
    return(res);
}

/* Jacobian matrix calculation */
jacobian(a, b, data) = {
    my(J, n, count, i);
    J = matrix(length(data), 2); \\ Initialize Jacobian matrix
    for (i = 1, length(data),
        n = data[i][2];
        count = data[i][1];
        J[i, 1] = -n^b; \\ Derivative with respect to a
        J[i, 2] = -a * n^b * log(n); \\ Derivative with respect to b
    );
    return(J);
}

/* Gauss-Newton method for fitting */
gauss_newton(data, a0, b0, max_iter, tol) = {
    my(a, b, iter, delta, J, r);
    a = a0; \\ Initialize a
    b = b0; \\ Initialize b

    for (iter = 1, max_iter,
        print("at iteration: ", iter);
        r = vector(length(data), i, data[i][2] - a * data[i][1]^b)~; \\ 
Residuals
        J = jacobian(a, b, data); \\ Compute Jacobian

        \\ Solve the normal equations J' * J * [delta_a, delta_b] = -J' * r
        my(JtJ, Jtr, scaled);
        JtJ = mattranspose(J) * J; \\ Compute J' * J
        Jtr = mattranspose(J) * r; \\ Compute J' * r
        scaled = 100000000; \\ initial guess

        delta = scaled * matsolve(JtJ, -Jtr); \\ Solve for delta
        print("delta : ",delta);

        \\ Update parameters
        a += delta[1]; \\ Update a
        b += delta[2]; \\ Update b

        \\ Optionally print current estimates
        print("Iteration: ", iter, " a: ", a, " b: ", b);

        \\ Check for convergence (could set a tolerance value)
        if (sqrt(delta[1]^2 + delta[2]^2) < tol, break);
    );

    return([a, b]); \\ Return the optimized parameters
}

resolve(file,tol) = {
    my(data, sz, a0, b0, max_iter, params);
    data = readvec(file);
    sz = matsize(data)[2];
    print("Input size: ", sz);
    a0 = 5.9791; \\ Initial guess for a
    b0 = 1.9288; \\ Initial guess for b
    max_iter = 100; \\ Maximum number of iterations

    \\ Run the Gauss-Newton fitting
    params = gauss_newton(data, a0, b0, max_iter, tol);
    print("Optimized parameters: ", params);
}