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.

To view all translated materials including this page, select Country 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
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.'