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

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);
}