Main Content

雅可比乘法函数与线性最小二乘法

打开实时脚本

使用雅可比乘法函数，可以求解如下形式的最小二乘问题

minx12‖C⋅x-d‖22

例如 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],

并且约束为 -5≤xi≤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，只需使用 v 的移位计算 C*Y，然后使用 fliplr(v) 的移位计算 C' 乘以结果。

辅助函数 lsqcirculant3 是一个实现此过程的雅可比乘法函数；它出现在本示例的末尾。

本例末尾的 dolsqJac3 辅助函数设置了向量v 并使用 lsqcirculant3 雅可比乘法函数调用了求解器 lsqlin。

当 n = 3000 时，C 是一个包含 18,000,000 个元素的密集矩阵。在选定的x值下，确定 n = 3000 的 dolsqJac3 函数的结果，并显示输出结构。

[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

另请参阅

|

相关主题