Используя функцию умножения якобиана, можно решить задачу вида методом наименьших квадратов
таким образом lb ≤ x ≤ ub
для задач, где C очень велик, возможно, слишком велик, чтобы храниться. Для этого метода используйте 'trust-region-reflective'
алгоритм.
Например, рассмотрим задачу, где C является 2n-на-n матрицей, основанной на циркулянтной матрице. Строки C являются сдвигами вектора-строки v. Этот пример имеет вектор-строку v с элементами вида :
,
где элементы смещены циклически.
Этот пример методом наименьших квадратов рассматривает задачу, где
,
и ограничениями являются для .
Для достаточно больших , плотная матрица C не помещается в память компьютера ( слишком велик в одной протестированной системе).
Функция умножения якобиана имеет следующий синтаксис.
w = jmfcn(Jinfo,Y,flag)
Jinfo
- матрица того же размера, что и C, используемая в качестве предкондиционера. Если C слишком велик, чтобы помещаться в память, Jinfo
должна быть разреженной. Y
является вектором или матрицей с таким размером, чтобы C*Y
или C'*Y
работает как матричное умножение. flag
говорит jmfcn
какой продукт сформировать:
flag
> 0 <reservedrangesplaceholder0>
flag
< 0 <reservedrangesplaceholder0>
flag
= 0 <reservedrangesplaceholder0>
Поскольку 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)
.
Функция помощника lsqcirculant3
- функция умножения якобиана, которая реализует эту процедуру; он появляется в конце этого примера.
The dolsqJac3
функция helper в конце этого примера настраивает вектор v и вызывает решатель lsqlin
использование lsqcirculant3
Функция умножения якобиана.
Когда n
= 3000, C является плотной матрицей с 18 000 000 элементов. Определите результаты dolsqJac3
функция для 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.
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.'
Этот код создает lsqcirculant3
вспомогательная функция.
function w = lsqcirculant3(Jinfo,Y,flag,v) % This function computes the Jacobian multiply function % 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 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
Этот код создает dolsqJac3
вспомогательная функция.
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 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