Correct use of algebra_solver? (getting odd error)


#1

I have a problem that requires a root finding algorithm. Originally I hard coded a Newton’s method algorithm, but with the extra operations it may make the code slower. So, I tried using the “algebra_solver” function instead:

functions
{
    vector algebra_system(vector y, vector theta, real[] x_r, int[] x_i)
    {
        vector[1] ret;
        
        ret[1] = y[1] - theta[1]*sin(y[1]) - theta[2];
        
        return ret;
    }
    
    real BoundMeanAnomaly(real e, real dt)
    {
        vector[2] theta;
        vector[1] y_guess;
        vector[1] ret;
        real x_r[0];
        int x_i[0];
        
        y_guess[1] = dt + e*sin(dt) + e*e*sin(2.0*dt)/2.0;
        theta[1] = e;
        theta[2] = dt;
        
        ret = algebra_solver(algebra_system, y_guess, theta, x_r, x_i);
        
        return ret[1];
    }
}

Here “BoundMeanAnomaly” used to have the Newton’s method code that I am replacing with the algebra_solver. But, I am getting these syntax errors:

second argument to algebra_solver (initial guess) must be data only and not reference parameters
fourth argument to algebra_solver (real data) must be data only and not reference parameters

It seems that “y_guess” must be data only (I’m guessing in the C++ code it needs to of type “double” and not type “stan::var”?). But, it depends on the inputed parameter e and dt! Also, with “x_r”, I have no clue on how to pass in a dummy “data only” vector.

I’m at a lost on how to implement this function.


#2

The solution (hopefully) shouldn’t depend on the initial guess though. There’s something in the pipe to make it possible to specify the initial guess as a function of parameters, but for now could you just make a constant guess? Like just y_guess[1] = 0?


#3

I may be able to get away with something like that (although I’d like to use this approximation since it get so close to the actual value), but even when I code this:

y_guess[1] = 0.0;

I still get the same error. It seems to doesn’t like y_guess being declared as a vector.


#4

What’s the error you get? Is it a compiler-time error or runtime error?


#5

It’s a parser error, here’s the full output (sorry, I should have included this in the original post):

Failed to parse Stan model ‘anon_model_00339554466c553c1c72af0c1d40877a’. Error message:
SYNTAX ERROR, MESSAGE(S) FROM PARSER:

second argument to algebra_solver (initial guess) must be data only and not reference parameters
fourth argument to algebra_solver (real data) must be data only and not reference parameters
error in ‘unkown file name’ at line 25, column 70

23:         theta[2] = dt;
24:         
25:         ret = algebra_solver(algebra_system, y_guess, theta, x_r, x_i);
                                                                         ^
26:         

PARSER EXPECTED: “)”


#6

you are running in some of the oddities of data inside functions…

what you need to do:

  1. create an argument y_guess to the BoundMeanAnomaly function
  2. define inside the transformed data block the y_guess and pass that into the BoundMeanAnomaly function

The problem is that any variable declared inside a function is internally propagated to be what Stan refers to a parameter - even if you assign a fixed constant to it. Doing what I describe above circumvents this problem.


#7

Well, that seems to have worked – thx! Although, it’ll be nice if it were possible to declare the “y_guess” input as a parameter. Is there anyway to base y_guess on the inputed variables (not data only variables) in the current interface?


#8

No, and that should be possible. The point would be that there wouldn’t be any gradient information propagated backwar based on the guess. It’d be easy to template out in the implementation with value_of. I created a Stan issue and assigned it to @charlesm93.


#9

Hi! There has been a version of the solver which can take y_guess as a parameter for a few months now. The fix will be available with Stan 2.18. You can also try to use a development version.


#10

Nice! I love it when issues I open are already solved.