Documentation

## Jacobian Multiply Function with Linear Least Squares

You can solve a least-squares problem of the form

`$\underset{x}{\mathrm{min}}\frac{1}{2}{‖C\cdot x-d‖}_{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 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=\left[\begin{array}{ccccc}1& -1/2& 1/3& ...& -1/n\\ -1/n& 1& -1/2& ...& 1/\left(n-1\right)\\ 1/\left(n-1\right)& -1/n& 1& ...& -1/\left(n-2\right)\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ -1/2& 1/3& -1/4& ...& 1\\ 1& -1/2& 1/3& ...& -1/n\\ -1/n& 1& -1/2& ...& 1/\left(n-1\right)\\ 1/\left(n-1\right)& -1/n& 1& ...& -1/\left(n-2\right)\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ -1/2& 1/3& -1/4& ...& 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.'```