Форум разработчиков электроники ELECTRONIX.ru > Арифметика высокой точности в ALGLIB, как ее использовать?

MSP430F

Dec 10 2013, 09:06

Всем доброго времени суток. Надеюсь, что мой вопрос не останется без ответа.
Использую функцию rmatrixsolvels из пакета alglib. Подключил просто - закинул в папку с проектом все хедеры и сишники alglib. Из 13 файлов пришлось подключить 6 для вызова лишь одной функции. exe подрос заметно в размерах, время компиляции увеличилось заметно, ну да ладно, главное цель достигнута. Решается система линейных уравнений с прямоугольной матрицей - число уравнений превышает число неизвестных. Минимизируется невязка. Но чувствую, не хватает точности. В alglib есть возможность работы с арифметикой высокой точности http://alglib.sources.ru/equations/linear.php Но как ее подключить ? Что-то я не соображу... Так все запутано в этом alglib. Использую компилятор MinGW-W64. Буду признателен за любые подсказки.

thermit

Dec 10 2013, 10:26

Цитата

MSP430F:
Но чувствую, не хватает точности.

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

TSerg

Dec 10 2013, 11:49

Проверьте с использованием мат.пакетов Матлаб, Маткад..

MSP430F

Dec 10 2013, 12:32

Цитата(thermit @ Dec 10 2013, 14:26)

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

Я проверил максимальную ошибку в вычислении свободных членов по найденным решениям, использую функции rmatrixsolve (для квадратной матрицы) и rmatrixsolvels (также в режиме квадратной матрицы). Ошибки отличаются на 3-5 порядков (в первой от 10-10 до 10-7, во второй от 10-5 до 10-3), хотя данные идентичны. Мне же предстоит решать систему не с квадратной, а прямоугольной матрицей. Если с квадратной такие ошибки,
то что же ожидать от прямоугольной ?

Xenia

Dec 10 2013, 12:47

Цитата(MSP430F @ Dec 10 2013, 16:32)

Я проверил максимальную ошибку в вычислении свободных членов по найденным решениям, использую функции rmatrixsolve (для квадратной матрицы) и rmatrixsolvels (также в режиме квадратной матрицы). Ошибки отличаются на 3-5 порядков (в первой от 10-10 до 10-7, во второй от 10-5 до 10-3), хотя данные идентичны. Мне же предстоит решать систему не с квадратной, а прямоугольной матрицей. Если с квадратной такие ошибки,
то что же ожидать от прямоугольной ?

Если у вас на квадратной матрице rmatrixsolvels дает ошибку от 10-5 до 10-3, то при переходе на прямоугольную матрицу хуже уже не будет. Ухудшение на "3-5 порядков" - это плата за метод более общего типа (в нем арифметики больше), которую вы уже заплатили.

thermit

Dec 10 2013, 13:47

Цитата

MSP430F:
Я проверил максимальную ошибку в вычислении свободных членов по найденным решениям, использую функции rmatrixsolve (для квадратной матрицы) и rmatrixsolvels (также в режиме квадратной матрицы). Ошибки отличаются на 3-5 порядков (в первой от 10-10 до 10-7, во второй от 10-5 до 10-3), хотя данные идентичны. Мне же предстоит решать систему не с квадратной, а прямоугольной матрицей. Если с квадратной такие ошибки,
то что же ожидать от прямоугольной ?

Исходные данные можете предоставить?

MSP430F

Dec 10 2013, 14:33

Цитата(Xenia @ Dec 10 2013, 16:47)

Если у вас на квадратной матрице rmatrixsolvels дает ошибку от 10-5 до 10-3, то при переходе на прямоугольною матрицу хуже уже не будет. Ухудшение на "3-5 порядков" - это плата за метод более общего типа (в нем арифметики больше), которую вы уже заплатили.

Спасибо за ответ! Я чувствовал нечто подобное, а Вы подтвердили мои догадки. Да, похоже, что использовать арифметику высокой точности по-любому очень проблемно здесь. Но если бы кто, подсказал, как, то можно было и попробовать. Если бы можно было просто заменить double на long double... Но там иакой финт не пройдет, слишком все в alglib наворочено.

Цитата(thermit @ Dec 10 2013, 17:47)

Исходные данные можете предоставить?

Да. В файле матрица 21x21 элемент и массив B из 21-го элемента (правая часть). Максимальная ошибка для этих данных соответственно у меня получилась 1,1х10-9 и 1,9х10-4.

thermit

Dec 10 2013, 15:31

1. Плохо обусловленная задача. Число обусловленности матрицы ~10^8. При погрешности входных данных 10^-16 ошибка в результате <=~10^-8
2. Боюсь, алглиб вас не спасет. Вам нужно свести исходную матрицу к треугольной при помощи последовательности отражений (преобразование хаусхолдера) или плоских вращений (преобразование гивенса) и найти решение обратной подстановкой. На вашем примере гивенс дает максимальную погрешность 1.5*10^-9.

ps
Ну или ковыряться с расширенной точностью, если ошибка 10^-8 вас не устраивает.

Xenia

Dec 10 2013, 16:27

Цитата(thermit @ Dec 10 2013, 19:31)

Боюсь, алглиб вас не спасет. Вам нужно свести исходную матрицу к треугольной при помощи последовательности отражений (преобразование хаусхолдера) или плоских вращений (преобразование гивенса) и найти решение обратной подстановкой. На вашем примере гивенс дает максимальную погрешность 1.5*10^-9.

Alglib действительно не спасет, там алгоритм для этого дела гадкий

. За Хаусхолдера же бояться нечего, это ортогональное преобразование, и ошибок оно практически не вносит (лишь на уровне точности арифметики). А вот на стадии обратной подстановки при плохо обусловленной матрице неминуемо возникнут проблемы, т.к. обусловленность полученной треугольной матрицы останется такой же малой, как у исходной.

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

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

MSP430F

Dec 10 2013, 16:44

Цитата(thermit @ Dec 10 2013, 19:31)

Вам нужно свести исходную матрицу к треугольной при помощи последовательности отражений (преобразование хаусхолдера) или плоских вращений (преобразование гивенса) и найти решение обратной подстановкой.

И где же взять код этих адгоритмов для прямоугольной матрицы, не подскажите? Для квадратной алгоритм вращений у меня есть, он очень неплох. А вот для прямоугольной...

Цитата(Xenia @ Dec 10 2013, 20:27)

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

Но ведь Alglib так и считает...

thermit

Dec 10 2013, 19:03

Цитата

MSP430F:
И где же взять код этих адгоритмов для прямоугольной матрицы, не подскажите? Для квадратной алгоритм вращений у меня есть, он очень неплох. А вот для прямоугольной...

Можно легко модифицировать для прямоугольной. Если исходники есть.

ps

Могу нарисовать за умеренную плату, гы...

MSP430F

Dec 11 2013, 07:51

Цитата(thermit @ Dec 10 2013, 23:03)

Можно легко модифицировать для прямоугольной. Если исходники есть.

Вот гляньте этот файл. Эта функция дает ошибку 2.7*10-8 на тех же исходных данных. Она не оригинальна, чуть модифицирована отсюда Мной добавлено создание копии матрицы в теле функции, а создание массива результатов убрано. Также введена проверка деления на 0.

thermit

Dec 11 2013, 11:00

Честно говоря, лень ковыряться с чужими исходниками.
Должно быть как-то так:
Ошибка ~9*10^-10

main.c

Код

#include <stdio.h>
#include <math.h>

int SolveByGivensRotation(double *A,double *b,int R,int C,double *x);

double a[21][21]=
{
{1.000000,    1.229402,    1.229566,    1.230016,    1.229648,    1.229771,    1.229975,    1.229361,    1.229525,    1.229648,    1.229484,    1.229238,    1.229156,    1.228992,    1.228214,    1.227518,    1.227928,    1.227108,    1.227231,    1.226453,    1.226576},
{1.000000,    1.229220,    1.228551,    1.229852,    1.229774,    1.230227,    1.230539,    1.230034,    1.229774,    1.229279,    1.228967,    1.229019,    1.228978,    1.228677,    1.228413,    1.227573,    1.227919,    1.227521,    1.226982,    1.226752,    1.226759},
{1.000000,    1.229566,    1.229689,    1.229846,    1.230118,    1.229974,    1.229979,    1.229705,    1.229582,    1.229385,    1.229313,    1.229417,    1.228537,    1.228060,    1.228113,    1.227828,    1.227652,    1.227244,    1.227276,    1.226674,    1.226855},
{1.000000,    1.228701,    1.229879,    1.229808,    1.230219,    1.230501,    1.230044,    1.229808,    1.229361,    1.229078,    1.229125,    1.229088,    1.228816,    1.228577,    1.227816,    1.228129,    1.227769,    1.227281,    1.227072,    1.227079,    1.226659},
{1.000000,    1.229696,    1.229674,    1.229618,    1.229730,    1.229786,    1.229718,    1.229696,    1.229607,    1.229741,    1.229528,    1.229237,    1.229349,    1.228286,    1.227962,    1.227481,    1.227738,    1.227279,    1.226708,    1.226932,    1.227044},
{1.000000,    1.229095,    1.229253,    1.230450,    1.229750,    1.229885,    1.229389,    1.229389,    1.229479,    1.229614,    1.229547,    1.229140,    1.229253,    1.228824,    1.228576,    1.227943,    1.227966,    1.227153,    1.226860,    1.226815,    1.226431},
{1.000000,    1.228654,    1.228773,    1.230379,    1.230498,    1.229963,    1.229784,    1.229190,    1.230022,    1.229784,    1.229071,    1.228773,    1.228952,    1.228297,    1.228714,    1.228119,    1.227822,    1.226810,    1.227286,    1.226929,    1.226989},
{1.000000,    1.229305,    1.229601,    1.229552,    1.229848,    1.230292,    1.229453,    1.229552,    1.229700,    1.229798,    1.229650,    1.228713,    1.228368,    1.229108,    1.228516,    1.227480,    1.228072,    1.227185,    1.227037,    1.227234,    1.226346},
{1.000000,    1.228883,    1.228629,    1.230443,    1.230170,    1.230034,    1.230053,    1.230190,    1.229663,    1.228766,    1.228766,    1.228961,    1.229000,    1.228746,    1.228473,    1.227440,    1.227849,    1.228122,    1.227089,    1.226640,    1.226893},
{1.000000,    1.229543,    1.229630,    1.229356,    1.229499,    1.229683,    1.230001,    1.229833,    1.229839,    1.229518,    1.229374,    1.229143,    1.229059,    1.229181,    1.228151,    1.227593,    1.227655,    1.227322,    1.227116,    1.226638,    1.226676},
{1.000000,    1.229429,    1.229542,    1.230312,    1.229936,    1.230049,    1.230105,    1.229336,    1.229148,    1.229054,    1.228716,    1.228491,    1.229636,    1.228359,    1.228453,    1.227590,    1.228266,    1.228040,    1.227609,    1.226276,    1.226464},
{1.000000,    1.229692,    1.229447,    1.229575,    1.229739,    1.230022,    1.229872,    1.229878,    1.229592,    1.229464,    1.229258,    1.229183,    1.229292,    1.228375,    1.227878,    1.227933,    1.227636,    1.227453,    1.227027,    1.227061,    1.226433},
{1.000000,    1.229186,    1.229563,    1.229939,    1.230316,    1.229971,    1.229469,    1.229312,    1.229374,    1.229531,    1.229469,    1.229280,    1.229061,    1.228527,    1.228810,    1.227617,    1.227335,    1.226644,    1.226833,    1.227021,    1.227554},
{1.000000,    1.229820,    1.229803,    1.230004,    1.229990,    1.229898,    1.229716,    1.229279,    1.229198,    1.229403,    1.229509,    1.229152,    1.229086,    1.228881,    1.227502,    1.227404,    1.227758,    1.227441,    1.227067,    1.226886,    1.227012},
{1.000000,    1.229545,    1.229452,    1.229675,    1.229582,    1.229861,    1.229805,    1.229712,    1.229694,    1.229824,    1.229285,    1.229174,    1.229025,    1.228413,    1.228394,    1.228227,    1.227707,    1.227206,    1.227002,    1.226575,    1.226649},
{1.000000,    1.229498,    1.229330,    1.229305,    1.229378,    1.229955,    1.230075,    1.229882,    1.229642,    1.229522,    1.229618,    1.229354,    1.228993,    1.228176,    1.228945,    1.227743,    1.227310,    1.227671,    1.226949,    1.226829,    1.226637},
{1.000000,    1.229748,    1.229748,    1.229999,    1.229971,    1.229776,    1.229608,    1.229385,    1.229469,    1.229497,    1.229329,    1.229134,    1.228799,    1.228687,    1.228576,    1.227516,    1.227237,    1.227432,    1.227460,    1.227153,    1.226288},
{1.000000,    1.228414,    1.229260,    1.229719,    1.229719,    1.230036,    1.230000,    1.229754,    1.229542,    1.229260,    1.229366,    1.229401,    1.229190,    1.228943,    1.228520,    1.228379,    1.228238,    1.226899,    1.226547,    1.226793,    1.226828},
{1.000000,    1.229178,    1.228977,    1.229153,    1.230486,    1.229706,    1.229857,    1.229304,    1.229304,    1.229405,    1.229556,    1.229480,    1.229027,    1.229153,    1.228675,    1.228398,    1.227694,    1.227719,    1.226814,    1.226487,    1.226436},
{1.000000,    1.229249,    1.229674,    1.229702,    1.229759,    1.229957,    1.229759,    1.229702,    1.229787,    1.229702,    1.229475,    1.229192,    1.228399,    1.227946,    1.227776,    1.227890,    1.228314,    1.227691,    1.228059,    1.226360,    1.226417},
{1.000000,    1.228960,    1.229996,    1.230061,    1.230319,    1.230190,    1.230190,    1.229219,    1.229414,    1.228960,    1.228766,    1.228572,    1.229090,    1.229155,    1.227990,    1.227860,    1.227666,    1.227278,    1.226890,    1.227278,    1.226955}
};

double b[21]=
{
1.226985,    1.226294,    1.226775,    1.225659,    1.227156,    1.226476,    1.225978,    1.225508,    1.226484,    1.225971,    1.226426,    1.226622,    1.224824,    1.226897,    1.226687,    1.226372,    1.226260,    1.226441,    1.226009,    1.226643,    1.226760
};

double aorig[21][21],aa[21][21],x[21],bb[21];

int main()
{
    int i,j;
    double acc,maxerror;

    for(i=0;i<21;i++)
    {
        for(j=0;j<21;j++) aa[i][j]=aorig[i][j]=a[j][i];
        bb[i]=b[i];
    }

SolveByGivensRotation(&aa[0][0],bb,21,21,x);

maxerror=0.0;
for(i=0;i<21;i++)
{
     acc=b[i];
     for(j=0;j<21;j++) acc-=aorig[i][j]*x[j];
     if(fabs(acc)>maxerror) maxerror=fabs(acc);
}

printf("Max error  = %.15g\n",maxerror);

return 0;

}

givens.c

Код

#include <math.h>

#define SQR(x)     (x*x)

int Rotation(double *A,double *b,int R,int C,int ri,int ci)
{
    int i;
    double t1,t2,t3,c,s;

    t1=A[C*ci+ci];
    t2=A[C*ri+ci];

    t3=sqrt(SQR(t1)+SQR(t2));
    if(t3==0.0) return 1;
    c=t1/t3;
    s=t2/t3;
    for(i=0;i<C;i++)
    {
        t1=c*A[ci*C+i]+s*A[ri*C+i];
        t2=c*A[ri*C+i]-s*A[ci*C+i];
        A[ci*C+i]=t1;
        A[ri*C+i]=t2;
    }
    t1=c*b[ci]+s*b[ri];
    t2=c*b[ri]-s*b[ci];
    b[ci]=t1;
    b[ri]=t2;

    return 0;
}

int GivensTriang(double *A,double *b,int R,int C)
{
    int i,j;

    for(i=0;i<C;i++)//For col
        for(j=R-1;j>i;j--) //For row
            if(Rotation(A,b,R,C,j,i)) return 1;

    return 0;
}

void BackwardSubstitution(double *A,double *b,double *x,int R,int C)
{
int i,j;
double tmp;

if(A[(C-1)*C+C-1]==0.0) x[C-1]=0.0;
else x[C-1]=b[C-1]/A[(C-1)*C+C-1];
for(i=C-2;i>=0;i--)
{
     tmp=b[i];
     for(j=i+1;j<C;j++) tmp-=x[j]*A[i*C+j];
     if(A[i*C+i]==0.0) x[i]=0.0;
     else x[i]=tmp/A[i*C+i];
}
}

int SolveByGivensRotation(double *A,double *b,int R,int C,double *x)
{
    if(GivensTriang(A,b,R,C)) return 1;
    BackwardSubstitution(A,b,x,R,C);

    return 0;
}

_pv

Dec 11 2013, 11:48

Код

void GaussSolve(double * A, double * Y, int num){
  //zeroeing left lower part of A[]
  for (int k = 0; k < num - 1; k++){   //diagonal
    for (int j = k + 1; j < num; j++){  //vertical
      double m = A[k + j * num] / A[k + k * num];
      for (int i = k+1; i < num; i++)  A[i + j * num] -= A[i + k * num] * m;  //horizontal
      Y[j] -= Y[k] * m;
    }
  }
  //zeroeing right upper part of A[], Y only.
  for (int k = num - 1; k > 0; k--)
    for (int j = k - 1; j >= 0; j--)
      Y[j] -= Y[k] * A[k + j * num] / A[k + k * num];
  //return answer in Y
  for (int i = 0; i < num; i++) Y[i] /= A[i + i * num];
}

обычный Гаусс на этих данных даёт ошибку 10^-12.
для прямоугольной матрицы при нахождении (A* A)^-1 A* поди не сильно больше наберётся.

thermit

Dec 11 2013, 12:07

Цитата

_pv:
обычный Гаусс на этих данных даёт ошибку 10^-12.

Не дает. Те же самые ~9*10^-10.
Матрицу нужно транспонировать.

MSP430F

Dec 11 2013, 13:01

Цитата(thermit @ Dec 11 2013, 15:00)

Честно говоря, лень ковыряться с чужими исходниками.
Должно быть как-то так:

Это Вы мне что, прямо таки код рабочий дали ?

Спасибо большое, буду проверять. Если все OK, то куда donation перевести ?

Так как же должно быть в BackwardSubstitution():

if(A[(C-1)*R+C-1]==0.0) x[C-1]=0.0;
else x[C-1]=b[C-1]/A[(C-1)*R+C-1]; или

if(A[(C-1)*C+C-1]==0.0) x[C-1]=0.0;
else x[C-1]=b[C-1]/A[(C-1)*C+C-1]; ?

_pv

Dec 11 2013, 13:29

5*10^-10
если то же самое в чём радость по сравнению с Гауссом?

thermit

Dec 11 2013, 14:20

Цитата(MSP430F @ Dec 11 2013, 16:01)

Это Вы мне что, прямо таки код рабочий дали ?

Спасибо большое, буду проверять. Если все OK, то куда donation перевести ?

Так как же должно быть в BackwardSubstitution():

if(A[(C-1)*R+C-1]==0.0) x[C-1]=0.0;
else x[C-1]=b[C-1]/A[(C-1)*R+C-1]; или

if(A[(C-1)*C+C-1]==0.0) x[C-1]=0.0;
else x[C-1]=b[C-1]/A[(C-1)*C+C-1]; ?

Последняя редакция верна. Первый раз скопировал с ошибкой.

Цитата

_pv:
если то же самое в чём радость по сравнению с Гауссом?

На плохо обусловленных системах с числом уравнений большим числа неизвестных ошибка гаусса будет на порядок больше. Нужно будет принимать специальные меры.
Кроме того, применение ортогонального преобразования гарантирует состоятельность и численную устойчивость алгоритма. Про гаусса такого сказать нельзя.
Если число операций не является приоритетом, лучше использовать вращения или отражения.

зы
старик уилкинсон утверждает, что ему всегда хватало гаусса. Но это было в эпоху медленных компьютеров...

TSerg

Dec 11 2013, 19:42

Цитата(thermit @ Dec 11 2013, 18:20)

старик уилкинсон утверждает, что ему всегда хватало гаусса. Но это было в эпоху медленных компьютеров...

Не только.
Все чаще дело имеем с разреженными системами, а Гаусс разреженность лучше сохраняет, в отличии от методов отражений и вращений.

MSP430F

Dec 12 2013, 09:01

Цитата(thermit @ Dec 11 2013, 18:20)

Я правильно догадался, что в этом месте ошибка:

for(j=0;j<21;j++) aa[i][j]=aorig[i][j]=a[j][i];

Должно быть:

for(j=0;j<21;j++) aa[i][j]=aorig[i][j]=a[i][j];

Первый индекс - строка, второй - столбец. Это же одинаково для всех массивов ?

thermit

Dec 12 2013, 09:18

Ошибки здесь нет. Исходная матрица лежит в памяти постолбцово (1-я строка состоит из единиц). Во всяком случае я так считал. Для работы алгоритмов матрица должна лежать в памяти построчно. for(j=0;j<21;j++) aa[i][j]=aorig[i][j]=a[j][i]; Это и есть такое преобразование (транспонирование).

MSP430F

Dec 12 2013, 10:18

Цитата(thermit @ Dec 12 2013, 13:18)

Ошибки здесь нет. Исходная матрица лежит в памяти постолбцово (1-я строка состоит из единиц). Во всяком случае я так считал. Для работы алгоритмов матрица должна лежать в памяти построчно. for(j=0;j<21;j++) aa[i][j]=aorig[i][j]=a[j][i]; Это и есть такое преобразование (транспонирование).

Если считать, что первый индекс - строка, а второй - столбец, то первая строка как раз и не состоит из 1.0, а равна
{1.000000, 1.229402, 1.229566, 1.230016, 1.229648, 1.229771, 1.229975, 1.229361, 1.229525, 1.229648, 1.229484, 1.229238, 1.229156, 1.228992, 1.228214, 1.227518, 1.227928, 1.227108, 1.227231, 1.226453, 1.226576},
как она и проинициализирована. Зачем же здесь тогда транспонирование ?

Можете проверить:

i = 0;
for(j=0; j<21; j++) printf("a = %f\n", a[i][j]);

_pv

Dec 12 2013, 10:28

тогда простой Гаусс даст ошибку в 10^-12.
а для прямоугольной матрицы тем же Гауссом постойте псевдообратную и будет вам счастье и без ALGLIBа.

thermit

Dec 12 2013, 10:56

Цитата

MSP430F:
Зачем же здесь тогда транспонирование ?

Для данных, которые вы прислали погрешность вычислялась (и с вашей совпадала) для матрицы, где строки являются столбцами из вашей таблицы. Поскольку я вставил ее себе в код копипастом, пришлось ввести транспонирование. С индексами все путем: 1-й - номер строки, 2-й - номер столбца.
Честно говоря, не понял проблемы. Разместите свои данные построчно да и все дела.

Цитата

_pv:
будет вам счастье и без ALGLIBа.

Собственно, то, что алглиб тут лишний, давно понятно.

MSP430F

Dec 12 2013, 12:03

Цитата(thermit @ Dec 12 2013, 14:56)

Для данных, которые вы прислали погрешность вычислялась (и с вашей совпадала) для матрицы, где строки являются столбцами из вашей таблицы. Поскольку я вставил ее себе в код копипастом, пришлось ввести транспонирование. С индексами все путем: 1-й - номер строки, 2-й - номер столбца.
Честно говоря, не понял проблемы. Разместите свои данные построчно да и все дела.

Нет, Вы не правы. В моих данных все на месте - строки это строки, столбцы - столбцы. Первый элемент каждой строки - это 1.0 (так и должно быть). Транспонирование не нужно.

thermit

Dec 12 2013, 12:15

Цитата

MSP430F:
Нет, Вы не правы. В моих данных все на месте - строки это строки, столбцы - столбцы. Первый элемент каждой строки - это 1.0 (так и должно быть). Транспонирование не нужно.

Это непринципиально. Кстати, что за данные? Откуда?

MSP430F

Dec 12 2013, 14:42

Цитата(thermit @ Dec 12 2013, 16:15)

Это непринципиально. Кстати, что за данные? Откуда?

Проверяю метод ближайших соседей для прогнозирования значения временного ряда хаотического типа. Для примера взял котировки EUR/USD с шагом 1 час.

Ваш код хорошо работает. Правда, пришлось немного попотеть, чтобы можно было в качестве аргумента передавать матрицу в виде нормального двумерного массива, так как размерность матрицы может меняться динамически. Ваш код дает ошибку не хуже, чем alglib. При небольшой прямоугольности разница заметна в пользу Вашего кода. Если матрица 21 * 25 или 21 * 31, то разница в ошибке незаметна. Зато код получается небольшой и нет кучи warning-ов при компиляции, как при использовании alglib.

thermit

Dec 12 2013, 14:49

Цитата

MSP430F:
Правда, пришлось немного попотеть, чтобы можно было в качестве аргумента передавать матрицу в виде нормального двумерного массива

Зачем передавать в виде двумерного массива? Решателю достаточно передать адрес начала матрицы и ее размерность.

Xenia

Dec 12 2013, 15:08

Цитата(thermit @ Dec 12 2013, 18:49)

Зачем передавать в виде двумерного массива? Решателю достаточно передать адрес начала матрицы и ее размерность.

У топик-стартера матрица аллокируется не одним блоком памяти, а построчно:

Код

  double **copy_mass = new double *[cnt_str]; // создаем копию матрицы, так как элементы матрицы будут модифицированы
  for(i = 0; i < cnt_str; i++){
    copy_mass[i] = new double [cnt_str + 1]; // N столбцов = N строк + 1
    for(j = 0; j < cnt_str + 1; j++) copy_mass[i][j] = mass[i][j];
  }

Однако это обстоятельство ничуть не мешает ее передавать через указатель double**, а еще лучше, как &double**.
Хотя на C++ можно было бы создать для матрицы класс, где кроме указателя double**, хранились бы и обе ее размерности. А конструктор и деструктор память под эту матрицу аллокировали и деаллокировали. В этом случае ссылку (здесь она удобнее указателя) на матрицу можно передавать одним единственным парамером. Кроме того, данный класс можно снабдить фукциями членами класса дял матричного умножения, сложения и т.п.