This is machine translation

Translated by Microsoft
Mouse over text to see original. Click the button below to return to the English verison of the page.

Jacobian Multiply Function with Linear Least Squares

You can solve a least-squares problem of the form


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:


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 shifted version of v times Y. The following matrix performs one step of the shift: v shifts to v*T, where


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)*T,fliplr(v)*T];

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 dolsqJac function in the following code sets up the vector v and matrix T, and calls the solver lsqlin:

function [x,resnorm,residual,exitflag,output] = dolsqJac(n)
r = 1:n-1; % index for making vectors

T = spalloc(n,n,n); % making a sparse circulant matrix
for m = r
T(n,1) = 1;

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 matrix T and vector v so they don't need to be
% computed in the Jacobian multiply function
options = optimoptions('lsqlin','Algorithm','trust-region-reflective',...

lb = -5*ones(1,n);
ub = 5*ones(1,n);

[x,resnorm,residual,exitflag,output] = ...

The Jacobian multiply function lsqcirculant is as follows:

function w = lsqcirculant(Jinfo,Y,flag,T,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);
    w = Jnegative(Jpositive(Y));

    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 = temp*T; % shift the circulant

    function a = Jnegative(q)
        % Calculate C'*q
        temp = fliplr(v)*T; % 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 = temp*T; % shift the circulant

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] = dolsqJac(3000);

Optimization terminated: relative function value changing by
 less than OPTIONS.FunctionTolerance.

ans =

ans =

ans =

output =
  struct with fields:
       iterations: 16
        algorithm: 'trust-region-reflective'
    firstorderopt: 7.5143e-05
     cgiterations: 36
          message: 'Optimization terminated: relative function value changing by les…'
Was this topic helpful?