Документация

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

Можно решить задачу наименьших квадратов формы

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

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

циклически переключенный:

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

и ограничениями является –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 который продукт сформироваться:

Поскольку 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.'

Смотрите также

circshift | fliplr

Похожие темы

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