Можно решить задачу наименьших квадратов формы
таким образом, что lb ≤ x ≤ ub, для проблем, где C является очень большим, возможно, слишком большим, чтобы храниться, при помощи якобиана, умножает функцию. Для этого метода используйте 'trust-region-reflective' алгоритм.
Например, рассмотрите случай, где C является матрицей 2n-by-n на основе циркулянтной матрицы. Это означает, что строки C являются сдвигами вектора-строки v. Этот пример имеет вектор-строку v с элементами формы (–1) k +1/k:
v = [1, –1/2, 1/3, –1/4..., –1/n],
циклически переключенный:
Этот пример наименьших квадратов рассматривает проблему где
d = [n – 1; n – 2;...; N,
и ограничениями является –5 ≤ x (i) ≤ 5 для i = 1..., n.
Для достаточно большого n плотный матричный C не помещается в память компьютера. (n = 10,000 является слишком большим в одной протестированной системе.)
Якобиан умножается, функция имеет следующий синтаксис:
w = jmfcn(Jinfo,Y,flag)
Jinfo матрица тот же размер как C, используемый в качестве предварительного формирователя. Если C является слишком большим, чтобы поместиться в память, Jinfo должно быть разреженным. Y вектор или матрица, измеренная так, чтобы C*Y или C'*Y целесообразен. flag говорит jmfcn какой продукт сформировать:
flag > 0 ⇒ w = C*Y
flag <0 ⇒ w = C'*Y
flag = 0 ⇒ w = C'*C*Y
Начиная с C такая просто структурированная матрица, легко записать, что якобиан умножает функцию в терминах векторного v; т.е. не формируя C. Каждая строка C*Y продукт циркулярной переключенной версии v времена Y. Используйте circshift к циркулярному сдвигу v.
Вычислить C*Y, вычислите v*Y чтобы найти первую строку, затем переключите v и вычислите вторую строку и так далее.
Вычислить C'*Y, выполните тот же расчет, но используйте переключенную версию temp, вектор, сформированный из первой строки C':
temp = [fliplr(v),fliplr(v)];
temp = [circshift(temp,1,2),circshift(temp,1,2)]; % Now temp = C'(1,:)Вычислить C'*C*Y, просто вычислите C*Y использование сдвигов v, и затем вычислите C' времена результат с помощью сдвигов fliplr(v).
dolsqJac3 функция в следующих кодовых наборах векторный v и вызывает решатель 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);Якобиан умножает функциональный lsqcirculant3 следующие:
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Когда n = 3000, C плотная матрица с 18,000,000 элементами. Вот результаты dolsqJac функция для n = 3000 в выбранных значениях x, и output структура:
[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.'