Re: fminunc and for loop

[Jordi Gutiérrez Hermoso]

On 10 July 2012 13:23, Michael Teitelbaum <address@hidden> wrote:
> I am trying to perform the fminunc function, but I want to use a for loop
> within my function, and I'm not sure exactly how to do that.
> What I want is:
>
> fminunc(f,[0,0])
> where f equals the sum of:
>
> f = @(x)((((A(i,3) / 1000 + x(1) + A(i,4) / 1000 - ((A(i,3) / 1000 + x(1) +
> A(i,4) / 1000) ^ 2 - 4 * A(i,3) / 1000 * A(i,4) / 1000) ^ 0.5) / 2 * 1.4e-3
> - (1.4e-3 - 8e-6) * (A(i-1,3) / 1000 + x(1) + A(i-1,4) / 1000 - ((A(i-1,3) /
> 1000 + x(1) + A(i-1,4) / 1000) ^ 2 - 4 * A(i-1,3) / 1000 * A(i-1,4) / 1000)
> ^ 0.5) / 2) * x(2)) * 1e6 - A(i,1)) ^ 2
>
> for i between 5 and  35.

Holy moly. Let's hope you don't have a typo in there...

Have you ever read this?

    
http://www.gnu.org/software/octave/doc/interpreter/Vectorization-and-Faster-Code-Execution.html

Use vector operations like vector indexing and sumsq. This is a
matrix-oriented language, after all.

At any rate, try something like the following:

    1;

    function out = midterm(A, idx, x)
      out =                                   \
          A(idx,3) + A(idx,4) + x(1)          \
          - ( (A(idx,3) + A(idx,4) + x(1)).^2 \
             - 4*A(idx,3).*A(idx,4))^.5;
    endfunction

    function out = f(A, x)

      ## Define indices of interest
      idx = 5:35;
      prev = idx-1;

      ## Scale the two columns that require it
      A(5:35,3:4) /= 1000;

      ## Compute the desired sum of squares objective function
      out = sumsq(1e6*x(2)*(midterm(A,idx)/2.4e-3
                            - 6.96e-4*midterm(A,prev))
                  - A(idx,1));
    endfunction

    result = fminunc(@(x) f(A,x), [0,0]);

Basically, write the common expression as another function
(copy-pasted code is a bug waiting to happen and a maintenance
nightmare), and use vector indexing and sumsq.

HTH,
- Jordi G. H.