Matrices in the objective function

[IgorTo]

The following code:

A = sprand(3000,3000,0.2);
x0 = rand(3000,1);

F = @(x) diag(x.^2)*A*x;
J_ad = AutoDiffJacobian(F, x0);

does not work. Error message is:

Requested 9000000x9000000 (73.0GB) array exceeds maximum
array size preference. Creation of arrays greater than this
limit may take a long time and cause MATLAB to become
unresponsive. See array size limit or preference panel for
more information.

Error in kron (line 37)
   K = matlab.internal.sparse.kronSparse(sparse(A),
   sparse(B));

Error in  *  (line 569)
                        My=kron(y',speye(size(x,1)));

It would be really good if this tool worked with objective functions which include manipulations with large sparse matrices.

[martinResearch]

For some reason matlab's kronSparse seems not to take into account the fact that the matrix will be sparse when computing the memory footprint of the resulting matrix. One option would be to rewrite the kronSparse function or a specialized function to compute

kron(a,speye(n))*b

however in your case it seems that simply adding parenteses to change the multiplication evaluation ordering solves the problem.
You can use

F = @(x) diag(x.^2)*(A*x);

or with a fix from the last commit you can avoid creating the matrix diag(x.^2) and use broadcasting (see here):

F = @(x) ((x.^2) .* A) * x;
F = @(x) (x.^2) .* (A * x);

The second line is faster to evaluate

[IgorTo]

Thank you so much! It is now very fast and does not return a memory error. Can't wait to try it in a Newton iteration for solving nonlinear elasticity PDEs.

[d-cogswell]

Hi, I'm running into a similar memory error while trying to solve a system of PDEs. The problem comes from the multiplication of NxN sparse matrices of the form:
A * spdiags(v, 0, N, N) * B

The call to spdiags() seems to be the issue. I replaced it with multiplications as you did above, but the performance is much worse. Any suggestions? Aside from the memory error with large simulations, I'm seeing good autodiff performance.

[martinResearch]

@d-cogswell , it would be great if you could provide me with some matrices A,B and vector v to replicate the problem. Did you try A*(v.*B)?

[d-cogswell]

@martinResearch sure, here's an example. I tried the alternate multiplication you suggested but it's really slow. Is it performing a dense multiply?

N=10000;
A=spdiags([-ones(N,1),ones(N,1)],0:1,N,N);
u=rand(N,1);

%This is fast, but fails with a memory error for large N
%N=20000 fails for me
f1 = @(x) A*spdiags(x,0,N,N)*A*x;
tic
J1=AutoDiffJacobianAutoDiff(f1, u);
toc

%This works for large N but it's really slow
f2 = @(x) A*(x.*A)*x;
tic
J2=AutoDiffJacobianAutoDiff(f2, u);
toc

[martinResearch]

Using a different ordering of the multiplications, it runs fast for me:

f3 = @(x) A*(x.*(A*x));
tic
J3=AutoDiffJacobianAutoDiff(f3, u);
toc

it runs in 0.16 sec on my machine for N=20000.

[d-cogswell]

Thanks for the tip! The order of operations really makes a big difference.