Ограничения алгоритма Левенберга-Марквардта

Question

Я использую алгоритм Левенберга-Марквардта , чтобы минимизировать нелинейную функцию 6 параметров.У меня есть около 50 точек данных для каждой минимизации, но я не получаю достаточно точных результатов.Может ли быть настолько значительным тот факт, что мои параметры отличаются друг от друга на несколько порядков?Если да, где мне искать решение?Если нет, какие ограничения LMA вы встретили в своей работе (это может помочь найти другие проблемы с моей заявкой)?Большое спасибо за вашу помощь.

Редактировать: Проблема, которую я пытаюсь решить, состоит в том, чтобы определить наилучшее преобразование T:

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

, чтобы установить набор трехмерных точек для группы трехмерных линий.Подробно у меня есть набор координат трехмерных точек и уравнения соответствующих трехмерных линий, которые должны проходить через эти точки (в идеальной ситуации).LMA минимизирует сумму расстояний преобразованных трехмерных точек до соответствующих трехмерных линий.Функция преобразования выглядит следующим образом:

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:

Некоторые примерные данные вставлены ниже.Трехмерные линии описываются центральной точкой и вектором направления.Центральная точка для всех линий равна (0,0,0), а координата 'uz' для каждого вектора равна 1. Набор координат 'ux' векторов направления:

-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

Набор 'uy'координаты векторов направления:

-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 ...):

 {{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

Что ж, правильный способ использования Левенберга-Марквардта состоит в том, что вам нужна хорошая начальная оценка («начальное число») для ваших параметров.Напомним, что LM является вариантом Ньютона-Рафсона;как и в случае с такими итеративными алгоритмами, качество вашей отправной точки сделает или сломает вашу итерацию;либо сходящиеся к тому, что вы хотите, сходящиеся к чему-то совершенно другому (что вряд ли произойдет, особенно если у вас много параметров), либо стрельба в дикий синий цвет (расходится).было бы более полезным, если бы вы могли упомянуть подходящую вам модель и, возможно, график разброса ваших данных;это может иметь большое значение в поиске подходящего решения для этого.

Answer 2

Здесь у вас есть проблема, смоделированная и работающая с Mathematica.

Я использовал метод Левенберга-Марквардта.

Вот почему я попросил ваши данные. С МОИМИ данными ВАШИ проблемы всегда будут легче:)

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

Из:

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

(в пределах 1/1000 от реальных значений)

Редактировать

Я немного поработал с твоими данными.
Проблема в том, что ваша система очень плохо подготовлена. Вам нужно гораздо больше данных, чтобы эффективно рассчитывать такие небольшие повороты.

Вот результаты, которые я получил:

Вращения в градусах:

rx = 179.99999999999999999999984968493536659553226696793
ry = 180.00000000000000000000006934755799995159952661222
rz = 180.0006286861217378980724139120849587855611645627

Traslations

tx = 48.503663696727576867196234527227830090575281353092
ty = 63.974139455057300403798198525151849767949596684232
tz = -0.99999999999999999999997957276716543927459921348549  

Я должен посчитать ошибки, но сейчас у меня нет времени.

Кстати, rz = Pi + 0,000011 (в радианах)

НТН!

Answer 3

Я бы посоветовал вам попробовать использовать другой подход для косвенного поиска параметров вращения, а именно использовать 4x4 матрицу аффинного преобразования для включения параметров перевода и вращения.

Это избавляетнелинейности функций синуса и косинуса (которые вы можете выяснить после факта).

Сложной частью будет ограничение матрицы преобразования от сдвига или масштабирования, что вам не нужно.

Answer 4

Ну, я использовал ceres-solver , чтобы решить эту проблему, но я внес изменения в ваши данные. Вместо «uz = 1.0» я использовал «uz = 0.0», что делает это полностью двумерным подбором данных.

Я получил следующие результаты. транс: -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];
}