Main Content

Using a Jacobian multiply function, 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. For this technique, use the `'trust-region-reflective'`

algorithm.

For example, consider a problem where *C* is a 2n-by-n matrix based on a circulant matrix. The rows of *C* are shifts of a row vector *v*. This example has the row vector *v* with elements of the form $$(\u20131{)}^{k+1}/k$$:

$$v=[1,-1/2,1/3,-1/4,\dots ,-1/n]$$,

where the elements are cyclically shifted.

$$C=\left[\begin{array}{ccccc}1& -1/2& 1/3& ...& -1/n\\ -1/n& 1& -1/2& ...& 1/(n-1)\\ 1/(n-1)& -1/n& 1& ...& -1/(n-2)\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -1/2& 1/3& -1/4& ...& 1\\ 1& -1/2& 1/3& ...& -1/n\\ -1/n& 1& -1/2& ...& 1/(n-1)\\ 1/(n-1)& -1/n& 1& ...& -1/(n-2)\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -1/2& 1/3& -1/4& ...& 1\end{array}\right].$$

This least-squares example considers the problem where

$$d=[n-1,n-2,\dots ,-n]$$,

and the constraints are $$-5\le {x}_{i}\le 5$$ for $$i=1,\dots ,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`

works as matrix multiplication. `flag`

tells `jmfcn`

which product to form:

`flag`

> 0 ⇒`w = C*Y`

`flag`

< 0 ⇒`w = C'*Y`

`flag`

= 0 ⇒`w = C'*C*Y`

Because *C* is such a simply structured matrix, you can easily write a Jacobian multiply function in terms of the vector *v*, 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 helper function `lsqcirculant3`

is a Jacobian multiply function that implements this procedure; it appears at the end of this example.

The `dolsqJac3`

helper function at the end of this example sets up the vector *v* and calls the solver `lsqlin`

using the `lsqcirculant3`

Jacobian multiply function.

When `n`

= 3000, *C* is an 18,000,000-element dense matrix. Determine the results of the `dolsqJac3`

function for `n`

= 3000 at selected values of *x*, and display 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.

disp(x(1))

5.0000

disp(x(1500))

-0.5201

disp(x(3000))

-5.0000

disp(output)

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

This code creates the `lsqcirculant3`

helper function.

function w = lsqcirculant3(Jinfo,Y,flag,v) % This function computes the Jacobian multiply function % 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 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

This code creates the `dolsqJac3`

helper function.

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 it 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 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); end