Существует ли библиотека для MQL4 по линейной алгебре? - страница 3

Новый комментарий
 
Mathemat wrote:
P.S. Кажется, я начинаю понимать, что вычисление определителя через миноры низших порядков - и правда задачка не по зубам чистому MQL4...


Дело не в МКЛ: в количестве и "вычислительной стоимости" операций - очень много операций деления для которых нужно держать максимально возможную точность представления чисел (в С это double, long double) иначе все заканчивается потерей точности, а не получением результата. На 1990 год вообще не существовало компьютеров способных в конечный период времени (менее тысячелетия, если не ошибаюсь) решать этим методом системы порядка выше 10 (цифры примерные - сейчас уже точно не помню, а лазить по старым лекциям лень). Вы это легко найдете в любом профессиональном курсе по численным методам.
В тоже время для реализации МКЭ (метод конечных элементов) мы ( кафедра Динамики и прочности машин. кафедра Прикладной математики Харьковского Политеха - тогда еще там работал) решали ситемы свыше 64К степеней свободы - там правда симметричные матрицы и разряженные (то есть не полностью заполненные), что несколько снижает требования к железу :).

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

Успехов.

ЗЫ Эксел решает методом Ньютона(если не ошибаюсь). Это итерационный метод, более дорогостоящий в плане вычислительных затрат, но более и универсальный - подходит не только для СЛАУ.
 
Спасибо, VladislavVG. У меня такой вопрос к Вам. Может ли при выборе метода решения как-то помочь тот факт, что все коэффициенты и свободные члены положительны и имеют небольшой разброс (ну, скажем, не более 15-20%)? Я понимаю, что есть вероятность слабой обусловленности матрицы. Вот потому и спрашиваю: Гаусс с выбором гл. элемента в этом случае адекватен - или лучше что-нибудь итерационное?

На 1990 год вообще не существовало компьютеров способных в конечный период времени (менее тысячелетия, если не ошибаюсь) решать этим методом системы порядка выше 10


Вы не ошиблись порядком? Уж больно маленькая циферка (правда, и год древний совсем для ЭВМ)...
 
Mathemat:
Спасибо, VladislavVG. У меня такой вопрос к Вам. Может ли при выборе метода решения как-то помочь тот факт, что все коэффициенты и свободные члены положительны и имеют небольшой разброс (ну, скажем, не более 15-20%)? Я понимаю, что есть вероятность слабой обусловленности матрицы. Вот потому и спрашиваю: Гаусс с выбором гл. элемента в этом случае адекватен - или лучше что-нибудь итерационное?

На 1990 год вообще не существовало компьютеров способных в конечный период времени (менее тысячелетия, если не ошибаюсь) решать этим методом системы порядка выше 10


Вы не ошиблись порядком? Уж больно маленькая циферка (правда, и год древний совсем для ЭВМ)...
Помочь может - вероятность работоспособности метода без выбора главного элемента (если брать Гаусса) возрастает. Но там, в методе с выбором, добавляется один иф и один цикл перестановок, что на скорость работы метода в целом (по сравнению с временем, необходимым для прямого и обратного хода) не существенно влияет. Но зато при совместных матрицах Вы 100% получите точное решение.
По поводу итерационных методов - не думаю, лучше возьмите Гаусса. Для того круга задач вполне и решение дает точное.
По поводу порядков - если и ошибся то не сильно. Кстати не факт, что сейчас машины позволяют решать намного большие задачи этим методом - там при возрастании порядка матрицы степенная (точнее, факториальная) зависимость количества операций. Условно: если тогда пороговая величина была 10, то сейчас 11-12, что, согласитесь, весьма немного.

Успехов.
 
Ну могу помочь с кодом для решения лин.системы:

#property copyright "Copyright v7 c 2006, Valery V. Chesnokov"
#property link      "http://chv.tele-kom.ru"
 
#include <stdlib.mqh>
 
extern int CountBars = 200;
 
// *** ИНДЕКСЫ в мартицах с 1-ЦЫ! ***
#define MatrixN 11
#define MatrixN1 12 // MUST: MatrixN + 1
#define MatrixN2 13 // MUST: MatrixN + 2
 
static double dMatrix[MatrixN1][MatrixN2]; // матрица (MatrixN; MatrixN + 1)
static double dMatrix_Sour[MatrixN1][MatrixN2]; // матрица (MatrixN; MatrixN + 1)
//int MatrixN = 11; // индексы: [0,..,MatrixN-1 ; 0,..,MatrixN]
static double dVectorSolve[MatrixN1]; // 1..MatrixN
 
static double Points;
 
static bool bSolved = false;
static bool bMatrixBusy = false;
 
// решить систему уравнений на матрице
bool SolveMatrix(int main_index)
{
   if (bMatrixBusy)
      return;
   bMatrixBusy = true;
 
   bool bRet = true;
   string fname;
   double d;
 
   // решить матрицу приведением к "треугольнику"
   for (int y = 1; y <= MatrixN; y++)
   {
      // одно преобразование к треугольному виду
      bRet = SolveOneCol(y);
      if (!bRet)
         break;
 
      d = y;
      fname = "matrix_s" + DoubleToStr(d, 0) + ".txt";
//      Print("Будет вызов WriteMatrixToFile, main_index = ", main_index, ", y = ", y, ", fname = ", fname);
      WriteMatrixToFile(fname, main_index, "Этапы решения");
   }   
 
   // получить единичную матрицу
   for (y = MatrixN; y >= 2; y--)
   {
      RowSolveUp(y);
 
      d = y;
      fname = "matrix_t" + DoubleToStr(d, 0) + "--.txt";
      WriteMatrixToFile(fname, main_index, "Треугольник вверх");
   }
 
   // взять решение
   for (int x = 1; x <= MatrixN; x++)
   {
      dVectorSolve[x] = dMatrix[x][MatrixN1];
   }
   
   CheckSolve();
 
   bMatrixBusy = false;
   
   return (bRet);
}
 
bool CheckSolve()
{
   // проверить решение //
   double d1, d2, diff, max_diff;
   
   max_diff = 0;
   
   for (int y = 1; y <= MatrixN; y++)
   {
      d1 = 0;
 
      for (int x = 1; x <= MatrixN; x++)
      {
         d1 = d1 + (dMatrix_Sour[y][x] * dVectorSolve[x]);
      }
      
      d2 = dMatrix_Sour[y][MatrixN1];
      diff = MathAbs(d1 - d2);
      Print("CheckSolve (y = ", y, "). d1 = ", d1, ", d2 = ", d2, ", diff = ", diff);
      
      if (max_diff < diff)
      {
         max_diff = diff;
      }
   }
 
   Print("CheckSolve max_diff = ", max_diff);
   
   return (true);
}
 
// данную строку умножить на свой коэф-нт, и провести её вверх
void RowSolveUp(int row)
{
   double d;
 
   if (row > 1)
   {
      for (int y = (row - 1); y >= 1; y--)
      {
         d = dMatrix[y][row];
 
         dMatrix[y][row] = 0;
         dMatrix[y][MatrixN1] = dMatrix[y][MatrixN1] - d * dMatrix[row][MatrixN1];
      }
   }
}
 
// одно преобразование матрицы к треугольному виду
bool SolveOneCol(int col)
{
   // найти max по модулю значение в столбце
   int row = col; // только нижняя треугольная часть матрицы
   double dMax = MathAbs(dMatrix[row][col]);
   
   for (int y = col; y <= MatrixN; y++)
   {
      if (dMax < MathAbs(dMatrix[y][col]))
      {
         dMax = MathAbs(dMatrix[y][col]);
         row = y;
      }
   }
   
   if (col != row)
      SwapRows(col, row, col);
 
   // разделить строку с max по модулю элементом на него
   if (!SolveDivideRow(col, col))
      return (false);
 
   // сделать во всех строках ниже в col столбце нули
   if (col < MatrixN)
   {
      for (y = col + 1; y <= MatrixN; y++)
      {
         SolveMultiPlusRow(col, col, y);
      }
   }
   
   return (true);
}
 
void SwapRows(int r1, int r2, int col)
{
   // поменять местами строки
   double d;
   
   for (int x = col; x <= MatrixN1; x++)
   {
      d = dMatrix[r1][x];
      dMatrix[r1][x] = dMatrix[r2][x];
      dMatrix[r2][x] = d;
   }
}
 
// разделить строку row, начиная с элемента col, на элемент matrix[row,col], он станет = 1.
bool SolveDivideRow(int row, int col)
{
   double d0 = dMatrix[row][col];
   
   if (d0 != 0)
   {
      dMatrix[row][col] = 1;
 
      // значения до col столбца уже равны 0
      for (int x = col + 1; x <= MatrixN1; x++)
      {
         dMatrix[row][x] = dMatrix[row][x] / d0;
      }
      
      return (true);
   }
   else
   {
      Print("SolveDivideRow(", row, ", ", col, ") дал 0 в ведущем элементе с max модулем. Матрица имеет бесконечное множество решений.");
      return (false);
   }
}
 
// прибавить к dest-строке другую sour-строку, предварительно умноженную н коэф-нт
void SolveMultiPlusRow(int col, int row_sour, int row_dest)
{
   double k = dMatrix[row_dest][col];
   
   if (k != 0)
   {
      dMatrix[row_dest][col] = 0;
      
      for (int x = col + 1; x <= MatrixN1; x++)
      {
         dMatrix[row_dest][x] = dMatrix[row_dest][x] - (dMatrix[row_sour][x] * k);
      }
   }
}
 
void WriteMatrixToFile(string filename, int main_index, string FirstMess)
{
//   if (bMatrixBusy)
//      return;
//   bMatrixBusy = true;
 
   int handle;
   datetime orderOpen = OrderOpenTime();
   handle = FileOpen(filename, FILE_CSV|FILE_WRITE, '\t');
   if (handle <= 0)
   {
      PrintErrMessage("FileOpen() FAILED ");
      bMatrixBusy = false;
      return;
   }
 
   string s;
 
   int i1 = main_index + 1;
   int i2 = main_index + MatrixN;
 
   FileWrite(handle, FirstMess, ": Valuta: " + Symbol() + ", Period: ", Period(), ", для предсказания " + main_index, "-го бара (", TimeToStr(Time[main_index]), ") были использованы бары от ", i1, " (", TimeToStr(Time[i1]), ") до ", i2, " (", TimeToStr(Time[i2]), ").");
 
   for (int y = 1; y <= MatrixN; y++)
   {
      s = "";
      for (int x = 1; x <= MatrixN1; x++)
      {
         s = s + DoubleToStr(dMatrix[y][x], 6) + " \t";
      }
      
      FileWrite(handle, s);
   }
 
   //FileWrite(handle, Close[0], Open[0], High[0], Low[0], TimeToStr(orderOpen));
   FileClose(handle);
 
//   bMatrixBusy = false;
}
 
// печатать сообщение об ошибке
// if (true), ERROR.
bool PrintErrMessage(string sStartMess)
{
   int err = GetLastError();
   if (err != 0)
   {
      // ErrorDescription возвращет описание
      Print(sStartMess, err, "; message: " , ErrorDescription(err));
      return (true);
   }
   
   return (false);
}
И вызов решения:
   WriteMatrixToFile("matrix_000_sour.txt", main_index, "Исходная матрица");
   bool bRet = SolveMatrix(main_index);
   WriteMatrixToFile("matrix_x_final.txt", main_index, "Преобразованная матрица");
   
   if (bRet == false)
      Print("Бесконечное множество решений.");
 
   Print("Матрица решена");
   bSolved = true;
 

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

P.S. Сам пытаюсь сотворить итерационный алгоритм по очень простой геометрической идее. Грубо говоря: выбираем произвольную точку в пространстве решений и строим от нее "перпендикуляры" ко всем гиперплоскостям, соответствующим уравнениям системы. Затем берем среднее арифметическое всех полученных векторов и получаем новую исходную точку, от которой снова строим "перпендикуляры". Оценка точности - по норме вектора ошибки между значением левой части системы в полученной точке и вектором свободного столбца. Так как над матрицей исходной системы мы при этом не производим никаких преобразований, то есть надежда, что процесс должен сходиться к решению независимо от степени обусловленности матрицы системы. Скорость сходимости, начиная с некоторого момента, похоже, стремится к экспоненциальной, хотя показатель экспоненты может не сильно отличаться от единицы, особенно если все гиперплоскости почти параллельны (т.е. система плохо обусловлена).

Получится - погоняю и сравню. У Вас, как я понял, Гаусс с выбором главного элемента?

 
Не совсем про линейные уравнения, но в ту же сторону. У интела есть пакет всевозможных математических подпрограмм (MKL). Скачать можно здесь:

http://www.intel.com/cd/software/products/asmo-na/eng/307757.htm

Я наткнулся на этот интеловский пакет после усердного поиска оптимизированных-по-скорости алгоритмов для разных математических проблем. Например, я нашёл следующее сравнение БПФ алгоритмов, в котором интеловский оказался победителем

http://www.fftw.org/speed/

Я скачал демо версию. Есть исходники кода, но всё по частям и разбросано.
 
Mathemat писал (а):


Получится - погоняю и сравню. У Вас, как я понял, Гаусс с выбором главного элемента?


Не могу 100% утверждать, не помню точно определение алгоритма Гаусса, но похоже, что так. Здесь матрица последовательно сводится к диагональной, и в обрабатываемом столбце выбирается элемент, максимальный по модулю. Как я помню из курса линейной алгебры, такой метод помогает в нахождении устойчивого решения. Хотя, как уже писал, после 12-15 порядка матрицы оно всё равно "расползалось", слишком велик был разброс порядков чисел в матрице после сведения к диагонали.
 

Драсьте Всем!

Я вот сколько не думал так и не смог придумать куда можно ввернуть метод решения системы линейных уровнения ?
Разве что только для метода наименьших квадратов только.
Просветите пожалуйста, а то уж больно интересно =)

 
chv wrote:
Хотя, как уже писал, после 12-15 порядка матрицы оно всё равно "расползалось", слишком велик был разброс порядков чисел в матрице после сведения к диагонали.

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

2 shobvas: например, в индикаторе полиномиальной регрессии в Code Base (см. код индикатора). Или в НС, состоящей из одного перцептрона с линейной функцией активации. Или в генетическом алгоритме с линейной аппроксимирующей функцией. Вообще говоря, к рынкету можно что угодно применить, совсем не обязательно зацикливаться на стандартных фильтрах и осцилляторах. Просто сам метод должен иметь статистическое основание для применения... А твой ник мне нравится :)
 
Ищю алгоритм вычисления обратной матрицы на MQL. Вот тут есть примеры как это делать http://alglib.sources.ru/matrixops/general/inv.php. Может кто то уже перевел на MQL. Заранее Спасибо.
1234
Новый комментарий
Причина обращения: