Кто-нибудь может потестировать функцию (линейная алгебра) ?

Candid 2008.10.10 12:17

Перевожу в MQL (попутно пытаясь оптимизировать) функцию LU декомпозиции матрицы и было бы замечательно, если бы кто-нибудь, располагающий работающим кодом с аналогичной своей функцией, попытался бы вставить вместо неё мою и сравнить результаты. Приз - моя функция(которая конечно моя только в части МКуЭЛизации) может оказаться быстрее и она не портит входную матрицу. То есть ориентирована на использование в индикаторах (можно на новом тике вычислять элементы А рекуррентно).

Описание:

A and n are input. LU is output, arranged as in equation

; indx[0..n-1] is an output vector that records the row permutation effected by the partial pivoting; d is output as +-1 depending on whether the number of row interchanges was even or odd, respectively.

A [0 ... nxn-1] - одномерный массив, содержит последовательно расположенные строки входной матрицы. LU - аналогичный выходной. То есть это квазидвумерные массивы.

Прежде чем использовать функцию, небходимо в подходящих местах сделать

  double TINY = 0.00001;
  TINY *= TINY;
  TINY *= TINY; // 1.0e-20, a small number.
  double vv[];           // vv stores the implicit scaling of each row.
  ArrayResize(vv,n);

И, наконец, сама функция

bool ludcmp(double A[], int n, double& LU[], int& indx[], double& d) {
  int i,imax,j,k;
  int ish,jsh,ksh,imaxsh;
  int i1,i2;
  double big,dum,sum,temp;
  int n1 = n-1;
  d=1.0;              // No row interchanges yet.
  ish = -n;
  for (i=0;i<n;i++) { // Loop over rows to get the implicit scaling information.
    ish += n;
    big=0.0;
    for (j=0;j<n;j++) {
      i2 = ish+j;
      temp = A[i2];
      LU[i2] = temp;
      if (temp < 0.0) temp = -temp;
      if (temp > big) big=temp;
    }
    if (big == 0.0) {
// No nonzero largest element.
      Print("function ludcmp: Singular matrix");
      return(false);
    } else {
      vv[i]=1.0/big;     // Save the scaling.
    }
  }
  jsh = -n;
  for (j=0;j<n;j++) {   // This is the loop over columns of Crout's method.
    jsh += n;
    ish = -n;
    for (i=0;i<j;i++) {  //  This is equation (2.3.12) except for i = j.
      ish += n;
      i2 = ish+j;
      sum=LU[i2];
      ksh = -n;
      for (k=0;k<i;k++) {
        ksh += n;
        sum -= LU[ish+k]*LU[ksh+j];
      }
      LU[i2]=sum;
    }
    big=0.0;             // Initialize for the search for largest pivot element.
    for (i=j;i<n;i++) {  // This is i = j of equation (2.3.12) and i = j+1 : ::N
      ish += n;
      i2 = ish+j;
      sum=LU[i2];          //  of equation (2.3.13). 
      ksh = -n;
      for (k=0;k<j;k++) {
        ksh += n;
        sum -= LU[ish+k]*LU[ksh+j];
      }
      LU[i2]=sum;
      if (sum < 0) sum = -sum;
      dum=vv[i]*sum;
      if (dum >= big) {
// Is the figure of merit for the pivot better than the best so far?
        big=dum;
        imax=i;
      }
    }
    if (j != imax) {     // Do we need to interchange rows?
      imaxsh = imax*n;
      for (k=0;k<n;k++) { // Yes, do so...
        i1 = imaxsh+k;
        i2 = jsh+k;
        dum=LU[i1];
        LU[i1]=LU[i2];
        LU[i2]=dum;
      }
      d = -d;          // ...and change the parity of d.
      vv[imax]=vv[j];      // Also interchange the scale factor.
    }
    indx[j]=imax;
    i2 = jsh+j;
    if (LU[i2] == 0.0) LU[i2]=TINY;
// If the pivot element is zero the matrix is singular (at least to the precision of the
// algorithm). For some applications on singular matrices, it is desirable to substitute
// TINY for zero.
    if (j != n1) {        // Now, finally, divide by the pivot element.
      dum=1.0/(LU[i2]);
      ish = i*j;
      for (i=j+1;i<n;i++) {
        ish += n;
        LU[ish+j] *= dum;
      }
    }
  }                    // Go back for the next column in the reduction.
  return(true);
}

Vladyslav Goshkov 2008.10.10 14:04 #1

lna01 >>:

А в чем проблема то ? Просто перемножьте матрицы - нижнюю (Л) на верхнюю (У), результат вычтите из исходной - все, что отличается от нуля больше чем на величину допустимой погрешности - ошибки.

Кстати, разложение Холецкого имеет и еще более привлекательную форму - LDU. Здесь в нижне- и верхне- трейгольных матрицах на диагоналях стоят 1, а D - диагональная матрица. В численных методах это существенно позволяет экономить временные затраты.

Успехов.

Вопрос математикам. Алгоритм вычисления Преобразования double Машинное обучение и нейронные

Candid 2008.10.10 15:05 #2

VladislavVG >>:

Прблемы нет, проверять я в любом случае буду. Но если найдётся более быстрая функция я предпочту использовать именно её :). Если конечно не удастся довести до нужных кондиций эту.

А разве разложение Холецкого не для частного случая симметричной матрицы?

Советники: Endless Possibilities Ошибки, баги, вопросы Любой вопрос новичка, чтоб

Vladyslav Goshkov 2008.10.10 16:19 #3

lna01 >>:

А разве разложение Холецкого не для частного случая симметричной матрицы?

Для несимметричных, если не запамятовал, называется обобщенное разложение Холецкого или LU (LDU как вариант).

Для Симметричной LDL(T). Или LL(T). Где L(T) - L транспонировання. Там просто считается одна матрица - вторая будет точно такой же. Реализован у меня такой алгоритм для профильных симметричных матриц. Еще в 1991 году делал на С для метода конечных элементов. Надо порыться в архивах. Сам алгоритм, если не ошибаюсь, представляет из себя прямой ход по Гауссу.

Успехов.

Существует ли библиотека для Некоторые признаки правильных ТС Создаем торгового робота

Candid 2008.10.10 17:28 #4

VladislavVG >>:

Для несимметричных, если не запамятовал, называется обобщенное разложение Холецкого или LU (LDU как вариант).

Успехов.

Сейчас у меня будут несимметрчные матрицы. Хотя идея где будут симметричные тоже есть, но она пока стоит в очереди и непонятно, когда она(очередь) подойдёт :). В принципе все нужные функции на Си у меня есть.

Кстати, Ваш прежний почтовый адрес (тот, что указывался в индикаторах) действителен? У меня есть один мелкий вопрос.

Есть ли ГРААЛЬ на Отображение объектов на всех Помогите разобраться в проблеме

Candid 2008.10.11 12:13 #5

Для проверки сделал на базе функции индикатор полиномиальной регрессии. Для сравнения использовал ang PR (Din)-v1. Для степеней, доступных ang PR (Din)-v1 результаты практически совпадают. Иногда наблюдаются пренебрежимые различия на хвостах, но это, думаю, объясяется деталями реализации. "Поднимаемая" степень зависит от размера выборки, для выборки в 720 точек максимальная степень - 53. Собственно для полинома это избыточно, но даёт представление о размерах систем, с которыми функция может работать.

Точное сравнение скоростей собственно процедур решения системы из-ра различий в алгоритмах было бы достаточно трудоёмким, грубая оценка показывает примерное равенство. Но, надо сказать, для регрессии матрица симметрична и для ускорения как раз можно пытаться использовать разложение Холецкого.

Таким образом считаю отладку функции законченной.

Файлы:

pr_lu.mq4 8 kb

ФР Н-волатильность Машинное обучение в трейдинге: Теория ускорения работы советника

Vladyslav Goshkov 2008.10.13 14:47 #6

lna01 >>:

Адрес действителен. 4vg AT mail.ru. Пишите, чем смогу помогу.

Успехов.

Candid 2008.10.13 16:55 #7

VladislavVG >>:

Адрес действителен. 4vg AT mail.ru. Пишите, чем смогу помогу.

Успехов.

Написал.

Vladyslav Goshkov 2008.10.13 17:06 #8

lna01 >>:

Написал.

Ответил. В смысле положительном. Пишу здесь на тот случай, если ответ потеряется.

Успехов.

Новый комментарий