This is machine translation

Translated by Microsoft
Mouseover text to see original. Click the button below to return to the English version of the page.

Note: This page has been translated by MathWorks. Click here to see
To view all translated materials including this page, select Country from the country navigator on the bottom of this page.

Jacobian Multiply Function with Linear Least Squares

You can solve a least-squares problem of the form

minx12Cxd22

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=[11/21/3...1/n1/n11/2...1/(n1)1/(n1)1/n1...1/(n2)1/21/31/4...111/21/3...1/n1/n11/2...1/(n1)1/(n1)1/n1...1/(n2)1/21/31/4...1].

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

See Also

|

Related Topics