2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки



Правила форума


В этом разделе нельзя создавать новые темы.

Если Вы хотите задать новый вопрос, то не дописывайте его в существующую тему, а создайте новую в корневом разделе "Помогите решить/разобраться (М)".

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

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

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



Начать новую тему Ответить на тему На страницу Пред.  1, 2, 3
 
 Re: МНК с особой точкой
Сообщение29.11.2014, 14:17 


28/08/14
9
Код для обращения таков (взят отсюда):
Код:
// Обращение квадратной матрицы
template<typename TT>
inline static bool InvertMN(int N, TT **M,
                            int         *GausOstatok=NULL,
                            long double *GausDet =   NULL,
                            long double *GausMinved= NULL) {
  // (c) Drobotenko http://www.drobotenko.com/bib_rus.html

  int gaus_ostatok;               // дефект линейного преобразования  0 == ОК
  long double gaus_deter;     // определитель в случае успешного обращения
  long double gaus_minved;  // минимальный ведущий элемент
                                           // можно использовать для оценки точности
  int *rn, *cn;
  int j, k;

  #define _NaN() (double&)*"Implementation dependent";
  //  заполнитель для неразрешенных измерений
  //  при невозможности обращения

  rn = new int [N];
  cn = new int [N];
  for(j=N; j--; )
    rn[j]=cn[j]=j;
  gaus_minved=1e99;
  gaus_deter=1;
  for(gaus_ostatok=N; gaus_ostatok; gaus_ostatok--) {
    int jved,kved;
    long double vved=-1, t;

    for(j=N; j--; ) {  // поиск ведущего
       if(~rn[j])
         for(k=N; k--; )
           if(~cn[k])
             if(vved<fabsl(M[j][k]))
                vved=fabsl(M[j][k]),jved=j,kved=k;  }

    if(gaus_minved>vved)
      gaus_minved=vved;
    gaus_deter*=M[jved][kved];
    if(vved<1e-99) {          // Для long double можно уменьшить
      for(j=N; j--; ) {
         if(~rn[j])
           for(k=N; k--; )
             M[j][k]=_NaN();  // ?
         if(~cn[j])
           for(k=N; k--; )
             M[k][j]=_NaN();  }
      delete[] rn;
      delete[] cn;
      return false;  }
    int jt=rn[jved],kt=cn[kved];
    // перестановки
    for(j=N; j--; )
      t=M[kt][j],M[kt][j]=M[jved][j],M[jved][j]=t;
    for(j=N; j--; )
      t=M[j][jt],M[j][jt]=M[j][kved],M[j][kved]=t;

    rn[jved]=rn[kt];
    cn[kved]=cn[jt];
    rn[kt]=cn[jt]=-1;

    vved=M[kt][jt];   M[kt][jt]=1;
    for(j=N; j--; )  {
      if(j==kt)
        continue;
      long double mul=M[j][jt]/vved;
      M[j][jt]=0;
      for(k=N; k--; )
        M[j][k]-=M[kt][k]*mul;  }  // самый внутренний цикл ровно N^3 операций
    for(k=N; k--; )
      M[kt][k]/=vved;
  }
  if (GausOstatok) *GausOstatok=gaus_ostatok;
  if (GausDet)        *GausDet =   gaus_deter;
  if (GausMinved)   *GausMinved= gaus_minved;

  delete[] rn;
  delete[] cn;
  return true;
}

Код работоспособный, при перемножении матриц $A^{-1}\cdot A$ размером 3x3 погрешности имеют порядок -19.
Матрицы, используемые в выражениях, заданы, предположительно, верно.
Выражение $(X^TX)^{-1} $ в случае $X=\begin{bmatrix}5 & 20 \\20 & 90 \end{bmatrix}$ вычисляется так:
$(X^TX)^{-1}=\begin{bmatrix}3,4 & -0,76 \\-0,76 & 0,17 \end{bmatrix}$,
аналогичные вычисления в MatLab дают те же результаты.

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 31 ]  На страницу Пред.  1, 2, 3

Модераторы: Модераторы Математики, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group