Функция умножения якобиана с линейными методами наименьших квадратов

Используя функцию умножения якобиана, можно решить задачу вида методом наименьших квадратов

minx12Cx-d22

таким образом     lb ≤ x ≤ ubдля задач, где C очень велик, возможно, слишком велик, чтобы храниться. Для этого метода используйте 'trust-region-reflective' алгоритм.

Например, рассмотрим задачу, где C является 2n-на-n матрицей, основанной на циркулянтной матрице. Строки C являются сдвигами вектора-строки v. Этот пример имеет вектор-строку v с элементами вида (1)k+1/k:

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

где элементы смещены циклически.

C=[1-1/21/3...-1/n-1/n1-1/2...1/(n-1)1/(n-1)-1/n1...-1/(n-2)-1/21/3-1/4...11-1/21/3...-1/n-1/n1-1/2...1/(n-1)1/(n-1)-1/n1...-1/(n-2)-1/21/3-1/4...1].

Этот пример методом наименьших квадратов рассматривает задачу, где

d=[n-1,n-2,,-n],

и ограничениями являются -5xi5 для i=1,,n.

Для достаточно больших n, плотная матрица C не помещается в память компьютера (n=10,000 слишком велик в одной протестированной системе).

Функция умножения якобиана имеет следующий синтаксис.

w = jmfcn(Jinfo,Y,flag)

Jinfo - матрица того же размера, что и C, используемая в качестве предкондиционера. Если C слишком велик, чтобы помещаться в память, Jinfo должна быть разреженной. Y является вектором или матрицей с таким размером, чтобы C*Y или C'*Y работает как матричное умножение. flag говорит jmfcn какой продукт сформировать:

  • flag0  <reservedrangesplaceholder0>

  • flag0  <reservedrangesplaceholder0>

  • flag0  <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

См. также

|

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте