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

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

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

Например, рассмотрим проблему, когда C является матрицей 2n-на-n на основе циркулирующей матрицы. Строки C являются сдвигами вектора строки v. Этот пример имеет вектор строки v с элементами вида $$(-1{)}^{k\mathrm{}}+$$1/k:

$$v=\mathrm{}[1\mathrm{,}\mathrm{-}\mathrm{}\mathrm{}\mathrm{1/2,1/3},\mathrm{}-1/4\mathrm{,}...$$, -1/n],

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

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

$$d=[\mathrm{}\mathrm{n-1},\mathrm{}n-2,..$$., -n],

и ограничения являются $$-\mathrm{}{}_{}\mathrm{}$$5≤xi≤5 для $$i\mathrm{}=1,$$..., n.

Для достаточно больших $$n$$плотная матрица C не помещается в память компьютера ($$n\mathrm{=}\mathrm{}\mathrm{}\mathrm{}\mathrm{}$$10000 слишком велика в одной тестируемой системе).

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

w = jmfcn(Jinfo,Y,flag)

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

Поскольку 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 - якобинская функция умножения, реализующая эту процедуру; в конце этого примера.

dolsqJac3 вспомогательная функция в конце этого примера устанавливает вектор v и вызывает решатель lsqlin с использованием lsqcirculant3 Якобская функция умножения.

Когда n = 3000, C - плотная матрица из 18 000 000 элементов. Определение результатов dolsqJac3 функция для n = 3000 при выбранных значениях x и отображение структуры вывода.

[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.'

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

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

См. также

circshift | fliplr

Связанные темы

• Квадратичная минимизация с плотным структурированным гессеном