Levenberg-Marquardtアルゴリズムの制限

Question

私はLevenberg-Marquardt algorithmを使用して、6つのパラメータの非線形関数を最小限に抑えています。私は各最小化のために約50データポイントを持っていますが、私は十分に正確な結果を得ていません。私のパラメータが数オーダーの大きさでお互いに異なるという事実は、それほど重要なことでしょうか？はいの場合は、どこで解決策を探すべきですか？そうでない場合は、あなたの仕事であなたが会ったLMAのどのような制限がありますか（それは私のapplictaionで他の問題を見つけるのに役立ちます）？ 多くのお手伝いをありがとうございます。

編集：3Dラインの束に3D点の集合にフィットするように

typedef struct 
{ 
    double x_translation, y_translation, z_translation; 
    double x_rotation, y_rotation, z_rotation; 
} transform_3D; 

：私が解決しようとしています問題は、最高の変換Tを決定することです。詳細には、3Dポイントの座標と対応する3Dラインの式のセットがあります。理想的な状況では、これらのポイントを通過する必要があります。 LMAは、変換された3D点の対応する3D線に対する距離の和を最小にする。 は、次のように関数である変換：...この説明は少し役立つ

cv::Point3d Geometry::transformation_3D(cv::Point3d point, transform_3D transformation) 
{ 
    cv::Point3d p_odd,p_even; 

    //rotation x 
    p_odd.x=point.x; 
    p_odd.y=point.y*cos(transformation.x_rotation)-point.z*sin(transformation.x_rotation); 
    p_odd.z=point.y*sin(transformation.x_rotation)+point.z*cos(transformation.x_rotation); 

    //rotation y 
    p_even.x=p_odd.z*sin(transformation.y_rotation)+p_odd.x*cos(transformation.y_rotation); 
    p_even.y=p_odd.y; 
    p_even.z=p_odd.z*cos(transformation.y_rotation)-p_odd.x*sin(transformation.y_rotation); 

    //rotation z 
    p_odd.x=p_even.x*cos(transformation.z_rotation)-p_even.y*sin(transformation.z_rotation); 
    p_odd.y=p_even.x*sin(transformation.z_rotation)+p_even.y*cos(transformation.z_rotation); 
    p_odd.z=p_even.z; 

    //translation 
    p_even.x=p_odd.x+transformation.x_translation; 
    p_even.y=p_odd.y+transformation.y_translation; 
    p_even.z=p_odd.z+transformation.z_translation; 

    return p_even; 
} 

希望を

をEDIT2：

いくつかの例示的なデータを以下に貼り付けられます。 3D線は、中心点および方向ベクトルによって記述される。全ての行のための中心点である（0,0,0）と「UZ」各ベクトルは方向ベクトルの「UX」座標の1 セットに等しい座標：「UY」座標の

-1.0986, -1.0986, -1.0986, 
-1.0986, -1.0990, -1.0986, 
-1.0986, -1.0986, -0.9995, 
-0.9996, -0.9996, -0.9995, 
-0.9995, -0.9995, -0.9996, 
-0.9003, -0.9003, -0.9004, 
-0.9003, -0.9003, -0.9003, 
-0.9003, -0.9003, -0.8011, 
-0.7020, -0.7019, -0.6028, 
-0.5035, -0.5037, -0.4045, 
-0.3052, -0.3053, -0.2062, 
-0.1069, -0.1069, -0.1075, 
-0.1070, -0.1070, -0.1069, 
-0.1069, -0.1070, -0.0079, 
-0.0079, -0.0079, -0.0078, 
-0.0078, -0.0079, -0.0079, 
0.0914, 0.0914, 0.0913, 
0.0913, 0.0914, 0.0915, 
0.0914, 0.0914 

セット

-0.2032, -0.0047, 0.1936, 
0.3919, 0.5901, 0.7885, 
0.9869, 1.1852, -0.1040, 
0.0944, 0.2927, 0.4911, 
0.6894, 0.8877, 1.0860, 
-0.2032, -0.0047, 0.1936, 
0.3919, 0.5902, 0.7885, 
0.9869, 1.1852, 1.0860, 
0.9869, 1.1852, 1.0861, 
0.9865, 1.1853, 1.0860, 
0.9870, 1.1852, 1.0861, 
-0.2032, -0.0047, 0.1937, 
0.3919, 0.5902, 0.7885, 
0.9869, 1.1852, -0.1039, 
0.0944, 0.2927, 0.4911, 
0.6894, 0.8877, 1.0860, 
-0.2032, -0.0047, 0.1935, 
0.3919, 0.5902, 0.7885, 
0.9869, 1.1852 

と（xyzxyzxyz ...）形式で3D点の集合：方向ベクトルの

{{0, 0, 0}, {0, 16, 0}, {0, 32, 0}, 
{0, 48, 0}, {0, 64, 0}, {0, 80, 0}, 
{0, 96, 0}, {0, 112,0}, {8, 8, 0}, 
{8, 24, 0}, {8, 40, 0}, {8, 56, 0}, 
{8, 72, 0}, {8, 88, 0}, {8, 104, 0}, 
{16, 0, 0}, {16, 16,0}, {16, 32, 0}, 
{16, 48, 0}, {16, 64, 0}, {16, 80, 0}, 
{16, 96, 0}, {16, 112, 0}, {24, 104, 0}, 
{32, 96, 0}, {32, 112, 0}, {40, 104, 0}, 
{48, 96, 0}, {48, 112, 0}, {56, 104, 0}, 
{64, 96, 0}, {64, 112, 0}, {72, 104, 0}, 
{80, 0, 0}, {80, 16, 0}, {80, 32, 0}, 
{80,48, 0}, {80, 64, 0}, {80, 80, 0}, 
{80, 96, 0}, {80, 112, 0}, {88, 8, 0}, 
{88, 24, 0}, {88, 40, 0}, {88, 56, 0}, 
{88, 72, 0}, {88, 88, 0}, {88, 104, 0}, 
{96, 0, 0}, {96, 16, 0}, {96, 32, 0}, 
{96, 48,0}, {96, 64, 0}, {96, 80, 0}, 
{96, 96, 0}, {96, 112, 0}} 

これは非常に小さな回転にデータをモデル化し、 "簡単" の一種です。

Answer 1

Levenberg-Marquardtを使用する適切な方法は、パラメータの初期推定値（シード値）が適切であることです。 LMはNewton-Raphsonの変形であることを思い出してください。このような反復アルゴリズムの場合と同様に、出発点の品質によって反復が行われます。あなたが望むものに収束するか、まったく異なるものに収束するか（特に、あなたが多くのパラメータを持っている場合は起こりそうもない）、または野生の青い雄鶏に発散する（発散する）。

あなたがフィッティングしているモデル関数と、おそらくデータの散布図について言えば、もっと役に立ちます。このためには実行可能な解決策を見出すために長い道のりを歩むかもしれません。

Answer 2

私は、回転パラメータを間接的に見つけるために別の方法、つまり4x4 affine transformation matrixを使用して変換パラメータと回転パラメータを組み込むことをお勧めします。

これは、正弦関数と余弦関数の非線形性を解消します（事実の後で分かります）。

難しい部分は、変換マトリックスが剪断またはスケーリングから拘束されることです。これは望ましくありません。

Answer 3

Mathematicaでモデル化して実行している問題があります。

私は "Levenberg-Marquardt"メソッドを使用しました。

私はあなたのデータを求めました。

{3.2423, {rx -> 0.500566, ry -> 0.334012, rz -> 0.199902, 
      tx -> 0.99985, ty -> 1.99939, tz -> 3.00021}} 

（実際の値の1000分の1以内）編集

：私のデータでは、あなたの問題は常に:)

xnew[x_, y_, z_] := 
    RotationMatrix[rx, {1, 0, 0}].RotationMatrix[ 
    ry, {0, 1, 0}].RotationMatrix[rz, {0, 0, 1}].{x, y, z} + {tx, ty, tz}; 

(* Generate Sample Data*) 
(* Angles 1/2,1/3,1/5 *) 
(* traslation -> {1,2,3} *) 
(* Minimum mean Noise 5% *) 

data = Table[{{x, y, z}, 
    RotationMatrix[1/2, {1, 0, 0}]. 
    RotationMatrix[1/3, {0, 1, 0}]. 
    RotationMatrix[1/5, {0, 0, 1}].{x, y, z} +{1, 2, 3} +RandomReal[{-.05, .05}, 3]}, 
    {x, 0, 1, .1}, {y, 0, 1, .1}, {z, 0, 1, .1}]; 

data = Flatten[data, 2]; 

(* Now find the parameters*) 
FindMinimum[ 
Sum[SquaredEuclideanDistance[xnew[i[[1]] /. List -> Sequence], 
    i[[2]]], {i, data}] 
, {rx, ry, rz, tx, ty, tz}, Method -> "LevenbergMarquardt"] 

アウトが容易になるだろう

私はあなたのデータで少し仕事をしました。
問題は、システムのコンディションが非常に悪いことです。そのような小さな回転を効果的に計算するには、さらに多くのデータが必要です。度で

ローテーション：

これらは私が得た結果である

rx = 179.99999999999999999999984968493536659553226696793 
ry = 180.00000000000000000000006934755799995159952661222 
rz = 180.0006286861217378980724139120849587855611645627 

Traslations

tx = 48.503663696727576867196234527227830090575281353092 
ty = 63.974139455057300403798198525151849767949596684232 
tz = -0.99999999999999999999997957276716543927459921348549 

私は、エラーを計算する必要がありますが、私は今は時間をしたん。

ところで、RZ =（ラジアン）パイ+ 0.000011

HTH！

Answer 4

私はこれを解決するためにceres-solverを使用しましたが、私はあなたのデータを変更しました。 "uz = 1.0"の代わりに、私は "uz = 0.0"を使用しました。これは完全に2dのデータフィッティングになります。

次の結果が得られました。 トランス：-88.6384、-16.3879、0 腐敗：0、0、-6.97813e-05

これらの結果を取得した後、手動で対応するラインに変換点の直交距離の和を計算し、0.0280452を得ました。

struct CostFunctor { 
    CostFunctor(const double p[3], double ux, double uy){ 
     p_[0] = p[0];p_[1] = p[1];p_[2] = p[2]; 
     n_[0] = ux; n_[1] = uy; 
     n_[2] = 0.0; 
     normalize(n_); 
    } 

    template <typename T> 
    bool operator()(const T* const x, T* residual) const { 
     T pDash[3]; 
     T pIn[3]; 
     T temp[3]; 
     pIn[0] = T(p_[0]); 
     pIn[1] = T(p_[1]); 
     pIn[2] = T(p_[2]); 
     //transform the input point p_ to pDash 
     xform(x, &pIn[0], &pDash[0]); 
     //find dot(pDash, n), where n is the direction of line 
     T pDashDotN = T(pDash[0]) * T(n_[0]) + T(pDash[1]) * T(n_[1]) + T(pDash[2]) * T(n_[2]); 
     //projection of pDash along line 
     temp[0] = pDashDotN * n_[0];temp[1] = pDashDotN * n_[1];temp[2] = pDashDotN * n_[2]; 
     //orthogonal vector from projection to point 
     temp[0] = pDash[0] - temp[0];temp[1] = pDash[1] - temp[1];temp[2] = pDash[2] - temp[2]; 
     //squared error 
     residual[0] = temp[0] * temp[0] + temp[1] * temp[1] + temp[2] * temp[2]; 
    return true; 
    } 
    //untransformed point 
    double p_[3]; 

    double ux_; 
    double uy_; 
    //direction of line 
    double n_[3]; 
}; 


template<typename T> 
void xform(const T *x, const T * inPoint, T *outPoint3) { 
    T xTheta = x[3]; 
    T pOdd[3], pEven[3]; 
    pOdd[0] = inPoint[0]; 
    pOdd[1] = inPoint[1] * cos(xTheta) + inPoint[2] * sin(xTheta); 
    pOdd[2] = -inPoint[1] * sin(xTheta) + inPoint[2] * cos(xTheta); 

    T yTheta = x[4]; 
    pEven[0] = pOdd[0] * cos(yTheta) + pOdd[2] * sin(yTheta); 
    pEven[1] = pOdd[1]; 
    pEven[2] = -pOdd[0] * sin(yTheta) + pOdd[2] * cos(yTheta); 


    T zTheta = x[5]; 

    pOdd[0] = pEven[0] * cos(zTheta) - pEven[1] * sin(zTheta); 
    pOdd[1] = pEven[0] * sin(zTheta) + pEven[1] * cos(zTheta); 
    pOdd[2] = pEven[2]; 

    T xTrans = x[0], yTrans = x[1], zTrans = x[2]; 
    pOdd[0] += xTrans; 
    pOdd[1] += yTrans; 
    pOdd[2] += zTrans; 

    outPoint3[0] = pOdd[0]; 
    outPoint3[1] = pOdd[1]; 
    outPoint3[2] = pOdd[2]; 
}