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 *A·x* ≤ *b*, *Aeq·x* = *beq*, *lb ≤ x ≤ ub*, for problems where *C* is
very large, perhaps too large to be stored, by using a Jacobian multiply
function.

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','JacobMult',... @(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.TolFun. x(1) ans = 5.0000 x(1500) ans = -0.5201 x(3000) ans = -5.0000 output output = iterations: 16 algorithm: 'trust-region-reflective' firstorderopt: 7.5143e-05 cgiterations: 36 message: 'Optimization terminated: relative function value changing by less...'

Was this topic helpful?