fminunc crashes matlab if large number of input variables used

I am working on a problem where input varibales with respect to which objective function (nonlinear) needs to be optimized is elements of a n*3 matrix where n is from 2000-3000. Here under I am providing essentional portion of the script.

clear all
clc
%n=2000;
n=3;k=3;normalized=0;periodic=0;
m=n+k+1;
[t,Range]=UniformKnotVector(k,n,periodic,normalized);

c=load('Data.txt');

p1=c(:,18:20);
p2=c(:,21:23);

v=zeros(size(c,1),3);

f=@(A)0;

for i=1:1:size(c,1) 

    disp(i)
    
    v(i,1)=(p2(i,1)-p1(i,1))/(sqrt(((p2(i,1)-p1(i,1)).^2)+((p2(i,2)-p1(i,2)).^2)+((p2(i,3)-p1(i,3)).^2)));
    v(i,2)=(p2(i,2)-p1(i,2))/(sqrt(((p2(i,1)-p1(i,1)).^2)+((p2(i,2)-p1(i,2)).^2)+((p2(i,3)-p1(i,3)).^2)));
    v(i,3)=(p2(i,3)-p1(i,3))/(sqrt(((p2(i,1)-p1(i,1)).^2)+((p2(i,2)-p1(i,2)).^2)+((p2(i,3)-p1(i,3)).^2)));
    
    x=c(i,24)*((n-k+2)/127.59);
    
    N11=BSpline1(k,t,x,n,periodic);
    N22=BSpline2(k,t,x,N11,periodic);
    N33=BSpline3(k,t,x,N22,periodic);
    
    fx=@(A)0;fy=@(A)0;fz=@(A)0;
    
    for j=1:1:(n+1)
        
       fx=@(A)(fx(A)+A(j,1)*N33(j,1));
       fy=@(A)(fy(A)+A(j,2)*N33(j,1));
       fz=@(A)(fz(A)+A(j,3)*N33(j,1));
       
    end
    
    f=@(A)(f(A)+((p1(i,1)-fx(A))^2+(p1(i,2)-fy(A))^2+(p1(i,3)-fz(A))^2-(((p1(i,1)-fx(A))*v(i,1))+((p1(i,2)-fy(A))*v(i,2))+((p1(i,3)-fz(A))*v(i,3)))^2));
    
end

A0=ones(n+1,3);

[Aopt,fopt]=fminunc(f,A0);

Script is working fine for smaller values of n (e.g., 3,4,5,...). With larger values (e.g., n=15,16,...) it is gradually taking more time. Finally for the actual n value for which I am interested (n=2000) script leads crash of MATLAB without any error.

What do I do?

Answer 1

[Matt J]

Well, the practice you have of cumulatively nesting anonymous functions inside one another is asking for trouble. You should at the very least convert to the whole thing to a non-anonymous function like below. Also, you are missing opportunities to do matrix/vector operations in place of for-loops.

   for i=1:1:size(c,1) 

    
    v(i,1)=(p2(i,1)-p1(i,1))/(sqrt(((p2(i,1)-p1(i,1)).^2)+((p2(i,2)-p1(i,2)).^2)+((p2(i,3)-p1(i,3)).^2)));
    v(i,2)=(p2(i,2)-p1(i,2))/(sqrt(((p2(i,1)-p1(i,1)).^2)+((p2(i,2)-p1(i,2)).^2)+((p2(i,3)-p1(i,3)).^2)));
    v(i,3)=(p2(i,3)-p1(i,3))/(sqrt(((p2(i,1)-p1(i,1)).^2)+((p2(i,2)-p1(i,2)).^2)+((p2(i,3)-p1(i,3)).^2)));
    
    x=c(i,24)*((n-k+2)/127.59);
    
    N11(:,i)=BSpline1(k,t,x,n,periodic);
    N22(:,i)=BSpline2(k,t,x,N11,periodic);
    N33(:,i)=BSpline3(k,t,x,N22,periodic);
    
   end
    
   
[Aopt,fopt]=fminunc(@(A)myObjective(A,N33,p,v),A0);

function f=myObjective(A,N33,p,v)  


    fxyz=N33.'*A;
    
    D=p1-fxyz;
    
    f=sum(D.^2-(D.*v).^2);

 end

[Matt J]

Note, your problem looks like an unconstrained linear least squares problem, and therefore likely has a closed-form, non-iterative solution.

Also, your objective looks separable in A(:,1), A(:,2), and A(:,3). It should be possible to solve for each separately.

[Satyajit Ghosh]

Dear Matt,

Thank you for quick answer. I have implemented your suggestions in a new script. It is running very fast now. But my optimized values of elements of matrix turns out to be erroneous.

Thanking you

Satyajit

[Satyajit Ghosh]

I checked the function value in two cases for a specific intial estimate matrix A0. The values are same in both cases.