複素数値データへのモデルの当てはめ

スクリプトを開く

この例では、複素数値データの非線形な当てはめを実行する方法を示します。ほとんどの Optimization Toolbox™ のソルバーとアルゴリズムは実数値データのみを処理しますが、最小二乗法ソルバーと fsolve は制約なしの問題の実数値データと複素数値データの両方を扱うことができます。目的関数は複素関数の観点から解析的でなければなりません。

複素数データを使用する場合には FunValCheck オプションを 'on' に設定しないでください。ソルバーでエラーが発生します。lsqcurvefit または lsqnonlin の 'interior-point' アルゴリズムを使用しないでください。このアルゴリズムは主に制約を処理するためのもので、複素数データでの動作は検証されていません。コード生成では、複素数データは処理されません。

データ モデル

データ モデルは単純な指数です。

$$y(x) = v_1 + v_2 e^{v_3 x}.$$

$x$ は入力データ、$y$ は応答、$v$ は係数の複素数値ベクトルです。ゴールは、$x$ とノイズを含む観測値 $y$ から $v$ を推定することです。データ モデルは解析的であるため、このモデルを複素数解で使用することができます。

ノイズを含む人為的なデータ

モデル用に人為的なデータを生成します。複素数係数ベクトル $v$ を [2;3+4i;-.5+.4i] とします。観測値 $x$ は指数分布しているとします。応答 $y$ に複素数値ノイズを追加します。

rng default % for reproducibility
N = 100; % number of observations
v0 = [2;3+4i;-.5+.4i]; % coefficient vector
xdata = -log(rand(N,1)); % exponentially distributed
noisedata = randn(N,1).*exp((1i*randn(N,1))); % complex noise
cplxydata = v0(1) + v0(2).*exp(v0(3)*xdata) + noisedata;

係数ベクトルを復元するためのモデルの当てはめ

データ モデルで予測された応答と観測値 ($x$ に対しては xdata、$y$ に対しては cplxydata) との差は次のとおりです。

objfcn = @(v)v(1)+v(2)*exp(v(3)*xdata) - cplxydata;

lsqnonlin または lsqcurvefit を使用してモデルをデータに当てはめます。この例では、まず lsqnonlin を使用します。

opts = optimoptions(@lsqnonlin,'Display','off');
x0 = (1+1i)*[1;1;1]; % arbitrary initial guess
[vestimated,resnorm,residuals,exitflag,output] = lsqnonlin(objfcn,x0,[],[],opts);
vestimated,resnorm,exitflag,output.firstorderopt

vestimated =

   2.1582 + 0.1351i
   2.7399 + 3.8012i
  -0.5338 + 0.4660i


resnorm =

  100.9933


exitflag =

     3


ans =

    0.0018

lsqnonlin によって、複素数係数ベクトルがほぼ有効数字 1 桁まで復元されます。残差のノルムが大きく、ノイズのためにモデルですべての観測値を当てはめられないことが示されています。1 次の最適性の尺度が約 1e-3 で、1e-6 未満ではないため、終了フラグは望ましい 1 ではなく 3 になっています。

代替方法: lsqcurvefit の使用

lsqcurvefit を使用して当てはめるには、応答から応答データを引いた値ではなく、応答のみが得られるようなモデルを作成します。

objfcn = @(v,xdata)v(1)+v(2)*exp(v(3)*xdata);

lsqcurvefit のオプションと構文を使用します。

opts = optimoptions(@lsqcurvefit,opts); % reuse the options
[vestimated,resnorm] = lsqcurvefit(objfcn,x0,xdata,cplxydata,[],[],opts)

vestimated =

   2.1582 + 0.1351i
   2.7399 + 3.8012i
  -0.5338 + 0.4660i


resnorm =

  100.9933

基礎となるアルゴリズムは同一であるため、結果は lsqnonlin と一致します。どちらでも、便利だと思うソルバーを使用してください。

代替方法: 実数部と虚数部の分割

範囲を含めたり、単に実数値の範囲内に完全に収めるには、係数の実数部と複素数部を別々の変数に分割できます。この問題では、係数を次のように分割します。

$$ \begin{array}{l}
y = {v_1} + i{v_2} + ({v_3} + i{v_4})\exp \left( {({v_5} + i{v_6})x} \right)\\
\ \ = \left( {{v_1} + {v_3}\exp ({v_5}x)\cos ({v_6}x) - {v_4}\exp ({v_5}x)\sin ({v_6}x)} \right)\\
\ \ + i \left( {{v_2} + {v_4}\exp ({v_5}x)\cos ({v_6}x) + {v_3}\exp ({v_5}x)\sin ({v_6}x)} \right).
\end{array}$$

lsqcurvefit の応答関数を作成します。

function yout = cplxreal(v,xdata)

yout = zeros(length(xdata),2); % allocate yout

expcoef = exp(v(5)*xdata(:)); % magnitude
coscoef = cos(v(6)*xdata(:)); % real cosine term
sincoef = sin(v(6)*xdata(:)); % imaginary sin term
yout(:,1) = v(1) + expcoef.*(v(3)*coscoef - v(4)*sincoef);
yout(:,2) = v(2) + expcoef.*(v(4)*coscoef + v(3)*sincoef);

このコードを MATLAB® パス上の cplxreal.m ファイルとして保存します。

応答データを実数部と虚数部に分割します。

ydata2 = [real(cplxydata),imag(cplxydata)];

係数ベクトル v が 6 つの次元をもつようになります。係数ベクトルをすべて 1 として初期化し、lsqcurvefit を使用して問題を解きます。

x0 = ones(6,1);
[vestimated,resnorm,residuals,exitflag,output] = ...
    lsqcurvefit(@cplxreal,x0,xdata,ydata2);
vestimated,resnorm,exitflag,output.firstorderopt

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.


vestimated =

    2.1582
    0.1351
    2.7399
    3.8012
   -0.5338
    0.4660


resnorm =

  100.9933


exitflag =

     3


ans =

    0.0018

6 要素のベクトル vestimated を 3 要素の複素数ベクトルとして解釈します。この解は以前の解と基本的に同じであることがわかります。

参考

トピック