You can solve a least-squares problem of the form
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 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:
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
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
endWhen 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 les…'