Jacobian Multiply Function with Linear Least Squares

You can solve a least-squares problem of the form

$$\underset{x}{\min }\frac{1}{2}{\Vert C\cdot x-d\Vert }_2^2$$

such that lb ≤ x ≤ ub, for problems where C is very large, perhaps too large to be stored, by using a Jacobian multiply function. For this technique, use the 'trust-region-reflective' algorithm.

For example, consider the case where C is a 2n-by-n matrix based on a circulant matrix. This means the rows of C are shifts of a row vector v. This example has the row vector v with elements of the form (–1)^k+1/k:

v = [1, –1/2, 1/3, –1/4, ... , –1/n],

cyclically shifted:

$$C=\left[\begin{array}{ccccc}1& -1/2& 1/3& \mathrm{...}& -1/n\\ -1/n& 1& -1/2& \mathrm{...}& 1/(n-1)\\ 1/(n-1)& -1/n& 1& \mathrm{...}& -1/(n-2)\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ -1/2& 1/3& -1/4& \mathrm{...}& 1\\ 1& -1/2& 1/3& \mathrm{...}& -1/n\\ -1/n& 1& -1/2& \mathrm{...}& 1/(n-1)\\ 1/(n-1)& -1/n& 1& \mathrm{...}& -1/(n-2)\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ -1/2& 1/3& -1/4& \mathrm{...}& 1\end{array}\right].$$

This least-squares example considers the problem where

d = [n – 1; n – 2; ...; –n],

and the constraints are –5 ≤ x(i) ≤ 5 for i = 1, ..., n.

For large enough n, the dense matrix C does not fit into computer memory. (n = 10,000 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 or C'*Y makes sense. flag tells jmfcn which product to form:

• flag > 0 ⇒  w = C*Y

• flag < 0 ⇒  w = C'*Y

• flag = 0 ⇒  w = C'*C*Y

Since C is such a simply structured matrix, it is easy to write a Jacobian multiply function in terms of the vector v; i.e., without forming C. Each row of C*Y is the product of a circularly shifted version of v times Y. Use circshift to circularly shift v.

To compute C*Y, compute v*Y to find the first row, then shift v and compute the second row, and so on.

To compute C'*Y, perform the same computation, but use a shifted version of temp, the vector formed from the first row of C':

temp = [fliplr(v),fliplr(v)];
temp = [circshift(temp,1,2),circshift(temp,1,2)]; % Now temp = C'(1,:)

To compute C'*C*Y, simply compute C*Y using shifts of v, and then compute C' times the result using shifts of fliplr(v).

The dolsqJac3 function in the following code sets up the vector v and calls the solver lsqlin:

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 that 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 really being 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);

The Jacobian multiply function lsqcirculant3 is as follows:

function w = lsqcirculant3(Jinfo,Y,flag,v)
% This function computes the Jacobian multiply functions
% 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% 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

When n = 3000, C is an 18,000,000-element dense matrix. Here are the results of the dolsqJac function for n = 3000 at selected values of x, and 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.

x(1)

ans =
    5.0000

x(1500)

ans =
   -0.5201

x(3000)

ans =
   -5.0000

output

output = 

  struct with fields:

       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.'