Documentation

### This is machine translation

Mouseover text to see original. Click the button below to return to the English version of the page.

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.

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

Watch now