Jacobian Multiply Function with Linear Least Squares

Open Live Script

Using a Jacobian multiply function, you can solve a least-squares problem of the form

$\min_{x} \frac{1}{2} {‖ C \cdot x - d ‖}_{2}^{2}$

such thatlb ≤ x ≤ ub, for problems whereCis very large, perhaps too large to be stored. For this technique, use the'trust-region-reflective'algorithm.

For example, consider a problem whereCis a 2n-by-n matrix based on a circulant matrix. The rows ofCare shifts of a row vectorv. This example has the row vectorvwith elements of the form $(- 1)^{k + 1} / k$ :

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

where the elements are cyclically shifted.

$C = [\begin{array}{cccccccccccccccccccc} 1 & - 1 / 2 & 1 / 3 & . . . & - 1 / n \\ - 1 / n & 1 & - 1 / 2 & . . . & 1 / (n - 1) \\ 1 / (n - 1) & - 1 / n & 1 & . . . & - 1 / (n - 2) \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ - 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) \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ - 1 / 2 & 1 / 3 & - 1 / 4 & . . . & 1 \end{array}] .$

This least-squares example considers the problem where

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

and the constraints are $- 5 \leq x_{i} \leq 5$ for $i = 1, \dots, n$ .

For large enough $n$ , the dense matrixCdoes 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)

Jinfois a matrix the same size asC, used as a preconditioner. IfCis too large to fit into memory,Jinfoshould be sparse.Yis a vector or matrix sized so thatC*YorC'*Yworks as matrix multiplication.flagtellsjmfcnwhich product to form:

flag> 0 ⇒w = C*Y
flag< 0 ⇒w = C'*Y
flag= 0 ⇒w = C'*C*Y

BecauseCis such a simply structured matrix, you can easily write a Jacobian multiply function in terms of the vectorv, without formingC. Each row ofC*Yis the product of a circularly shifted version ofvtimesY. Usecircshiftto circularly shiftv.

To computeC*Y, computev*Yto find the first row, then shiftvand compute the second row, and so on.

To computeC'*Y, perform the same computation, but use a shifted version oftemp, the vector formed from the first row ofC':

temp = [fliplr(v),fliplr(v)];

temp = [circshift(temp,1,2),circshift(temp,1,2)]; % Now temp = C'(1,:)

To computeC'*C*Y, simply computeC*Yusing shifts ofv, and then computeC'times the result using shifts offliplr(v).

The helper functionlsqcirculant3is a Jacobian multiply function that implements this procedure; it appears at theend of this example.

ThedolsqJac3helper function at theend of this example建立了向量vand calls the solverlsqlinusing thelsqcirculant3Jacobian multiply function.

Whenn= 3000,Cis an 18,000,000-element dense matrix. Determine the results of thedolsqJac3function forn= 3000 at selected values ofx, 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.'

Helper Functions

这段代码创建了lsqcirculant3helper function.

functionw = lsqcirculant3(Jinfo,Y,flag,v)% This function computes the Jacobian multiply function% for a 2n-by-n circulant matrix example.ifflag > 0 w = Jpositive(Y);elseif国旗< 0 w = Jnegative (Y);elsew = Jnegative(Jpositive(Y));endfunctiona = Jpositive(q)% Calculate C*qtemp = v; a = zeros(size(q));% Allocating the matrix aa = [a;a];% The result is twice as tall as the input.forr = 1:size(a,1) a(r,:) = temp*q;% Compute the rth rowtemp = circshift(temp,1,2);% Shift the circulantendendfunctiona = Jnegative(q)% Calculate C'*qtemp = 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 aforr = 1:len a(r,:) = [temp,temp]*q;% Compute the rth rowtemp = circshift(temp,1,2);% Shift the circulantendendend

这段代码创建了dolsqJac3helper function.

function[x,resnorm,residual,exitflag,output] = dolsqJac3(n)%r = 1:n-1;% Index for making vectorsv(n) = (-1)^(n+1)/n;% Allocating the vector vv(r) =( -1).^(r+1)./r;%现在应该2 n - Cby-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

Jacobian Multiply Function with Linear Least Squares

Helper Functions

See Also

Related Topics