Jacobian Multiply Function with Linear Least Squares
Using a Jacobian multiply function, you can solve a least-squares problem of the form
lb ≤ x ≤ ub, for problems where C is very large, perhaps too large to be stored. For this technique, use the
For example, consider a problem where C is a 2n-by-n matrix based on a circulant matrix. The rows of C are shifts of a row vector v. This example has the row vector v with elements of the form :
where the elements are cyclically shifted.
This least-squares example considers the problem where
and the constraints are for .
For large enough , the dense matrix C does not fit into computer memory ( is too large on one tested system).
A Jacobian multiply function has the following syntax.
w = jmfcn(Jinfo,Y,flag)
Jinfo is a matrix the same size as C, used as a preconditioner. If C is too large to fit into memory,
Jinfo should be sparse.
Y is a vector or matrix sized so that
C'*Y works as matrix multiplication.
jmfcn which product to form:
flag> 0 ⇒
w = C*Y
flag< 0 ⇒
w = C'*Y
flag= 0 ⇒
w = C'*C*Y
Because C is such a simply structured matrix, you can easily write a Jacobian multiply function in terms of the vector v, without forming C. Each row of
C*Y is the product of a circularly shifted version of v times
circshift to circularly shift v.
v*Y to find the first row, then shift v and compute the second row, and so on.
C'*Y, perform the same computation, but use a shifted version of
temp, the vector formed from the first row of
temp = [fliplr(v),fliplr(v)];
temp = [circshift(temp,1,2),circshift(temp,1,2)]; % Now temp = C'(1,:)
C'*C*Y, simply compute
C*Y using shifts of v, and then compute
C' times the result using shifts of
The helper function
lsqcirculant3 is a Jacobian multiply function that implements this procedure; it appears at the end of this example.
dolsqJac3 helper function at the end of this example sets up the vector v and calls the solver
lsqlin using the
lsqcirculant3 Jacobian multiply function.
n = 3000, C is an 18,000,000-element dense matrix. Determine the results of the
dolsqJac3 function for
n = 3000 at selected values of x, and display the output structure.
[x,resnorm,residual,exitflag,output] = dolsqJac3(3000);
Local minimum possible. lsqlin stopped because the relative change in function value is less than the function tolerance.
iterations: 16 algorithm: 'trust-region-reflective' firstorderopt: 5.9351e-05 cgiterations: 36 constrviolation:  linearsolver:  message: 'Local minimum possible.↵↵lsqlin stopped because the relative change in function value is less than the function tolerance.'
This code creates the
lsqcirculant3 helper function.
function w = lsqcirculant3(Jinfo,Y,flag,v) % This function computes the Jacobian multiply function % for a 2n-by-n circulant matrix example. if flag > 0 w = Jpositive(Y); elseif flag < 0 w = Jnegative(Y); else w = Jnegative(Jpositive(Y)); end function a = Jpositive(q) % Calculate C*q temp = v; a = zeros(size(q)); % Allocating the matrix a a = [a;a]; % The result is twice as tall as the input. for r = 1:size(a,1) a(r,:) = temp*q; % Compute the rth row temp = circshift(temp,1,2); % Shift the circulant end end function a = Jnegative(q) % Calculate C'*q temp = fliplr(v); temp = circshift(temp,1,2); % Shift the circulant for C' len = size(q,1)/2; % The returned vector is half as long % as the input vector. a = zeros(len,size(q,2)); % Allocating the matrix a for r = 1:len a(r,:) = [temp,temp]*q; % Compute the rth row temp = circshift(temp,1,2); % Shift the circulant end end end
This code creates the
dolsqJac3 helper function.
function [x,resnorm,residual,exitflag,output] = dolsqJac3(n) % r = 1:n-1; % Index for making vectors v(n) = (-1)^(n+1)/n; % Allocating the vector v v(r) =( -1).^(r+1)./r; % Now C should be a 2n-by-n circulant matrix based on v, % but it might be too large to fit into memory. r = 1:2*n; d(r) = n-r; Jinfo = [speye(n);speye(n)]; % Sparse matrix for preconditioning % This matrix is a required input for the solver; % preconditioning is not used in this example. % Pass the vector v so that it does not need to be % computed in the Jacobian multiply function. options = optimoptions('lsqlin','Algorithm','trust-region-reflective',... 'JacobianMultiplyFcn',@(Jinfo,Y,flag)lsqcirculant3(Jinfo,Y,flag,v)); lb = -5*ones(1,n); ub = 5*ones(1,n); [x,resnorm,residual,exitflag,output] = ... lsqlin(Jinfo,d,,,,,lb,ub,,options); end