Note: This page has been translated by MathWorks. Please click here

To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

You can solve a least-squares problem of the form

$$\underset{x}{\mathrm{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 2*n*-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)\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -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)\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -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 shifted version of `v`

times `Y`

.
The following matrix performs one step of the shift: `v`

shifts
to `v*T`

, where

$$T=\left[\begin{array}{ccccc}0& 1& 0& \mathrm{...}& 0\\ 0& 0& 1& \mathrm{...}& 0\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0& 0& 0& \mathrm{...}& 1\\ 1& 0& 0& \mathrm{...}& 0\end{array}\right].$$

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(m,m+1)=1; end 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',... 'JacobianMultiplyFcn',@(Jinfo,Y,flag)lsqcirculant(Jinfo,Y,flag,T,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 `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); 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 = temp*T; % shift the circulant end end 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 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] = dolsqJac(3000); Optimization terminated: relative function value changing by less than OPTIONS.FunctionTolerance. 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: 7.5143e-05 cgiterations: 36 message: 'Optimization terminated: relative function value changing by les…'

Was this topic helpful?