Main Content

線形最小二乗付きヤコビ乗算関数

ヤコビ乗算関数を使用すれば、次の形式の最小二乗問題を解くことができます。

minx12Cx-d22

これは lb ≤ x ≤ ub などが条件になり、"C" が非常に大きく格納できない問題です。この手法に対しては、'trust-region-reflective' アルゴリズムを使用してください。

たとえば、C が循環行列をベースにした 2n 行 n 列の行列である問題を考えます。C の行は行ベクトル v をシフトしたものです。この例は、(1)k+1/k 形式の要素をもつ次の行ベクトル v をもちます。

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

制約は i=1,,n に対して -5xi5 です。

n が非常に大きくなると、密行列 C はコンピューターのメモリに収まりません (n=10,000 は 1 つのテスト システムに対しては大きすぎます)。

ヤコビ乗算関数の構文は以下になります。

w = jmfcn(Jinfo,Y,flag)

JinfoC と同じサイズの行列であり、前提条件子として使用されます。C がメモリに収めるには大きすぎる場合、Jinfo をスパース行列にします。YC*Y または C'*Y が行列乗算として機能するようにサイズ設定されたベクトルまたは行列です。flagjmfcn に以下のような乗算の形式を指定します。

  • flag > 0 ⇒ w = C*Y

  • flag < 0 ⇒ w = C'*Y

  • flag = 0 ⇒ w = C'*C*Y

C はシンプルな構造の行列であるため、ベクトル v のヤコビ乗算関数は、C の作成を必要とせず容易に記述できます。C*Y の各行は vY の積を循環的にシフトしたものの積です。circshift を使用して、v を循環的にシフトします。

C*Y を計算するには、1 番目の行を見つけるために v*Y を計算し、次に v をシフトして 2 番目の行を計算する、というようにします。

C'*Y の計算には同じ計算を実行しますが、temp をシフトしたものを使用します。これは C' の 1 番目の行から形成されます。

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 (この例の終わりに掲載) は、ヤコビ乗算関数 lsqcirculant3 を使用して、ベクトル v を設定し、ソルバー 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

参考

|

関連するトピック