[slogos]

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

[kp+]

Чувствую, что это как то связано с числом обусловленности, надо бы исследовать эту тему. Подскажите, пж-ста где почерпнуть доступную для не математика информацию.

Жжош! Нематематики типа меня даже слов таких не знают
А если серьезно, в центре ресурсов (Resource Center) MathCAD'а много интересного, причем готового к употреблению. И в пакете Mathemathica от Wolfram Research тоже есть подобный раздел. Может, там все и найдется

[Сергей Дубина]

Если ответ нашел, ответь плз.
Одно время тоже бился над этим вопросом (но без этих страшных терминов) вот образец попытки в эксэль (матрица, транспонированная, произведение, обратная (ошибка))
-2 -6
0 3
1 8

-2 0 1
-6 3 8

40 -18 -50
-18 9 24
-50 24 65

-3,97E+13#ЗНАЧ! #ЗНАЧ!
#ЗНАЧ! #ЗНАЧ! #ЗНАЧ!
#ЗНАЧ! #ЗНАЧ! #ЗНАЧ!

З.Ы. Свою задачу я бросил, и решил ее по другому. Если кто даст разъяснения, буду говорить спасибо.

[slogos]

Прорешал пример, обратная получилась, но коэффициенты аж в 14 степени!
Еденичная, естесственно "плохая" - диагональные элементы далеко не равны 1.
Вопрос и заключается в том, чем объяснить, что при "вроде бы нормальных" коэффициентах исходной матрицы коэффициенты обратной такие "невообразимые".
Мне не удалось найти источник с объяснением этого факта.
И второе как из такого "выгребать", ведь если коэффициенты недостоверные, то и искомые величины будут недостоверными.

[kp+]

http://forum.exponenta.ru/ - здесь можно спросить у более близких к математике людей

[Елпанов Евгений]

Сообщение от slogos Прорешал пример, обратная получилась, но коэффициенты аж в 14 степени! Еденичная, естесственно "плохая" - диагональные элементы далеко не равны 1. Вопрос и заключается в том, чем объяснить, что при "вроде бы нормальных" коэффициентах исходной матрицы коэффициенты обратной такие "невообразимые". Мне не удалось найти источник с объяснением этого факта. И второе как из такого "выгребать", ведь если коэффициенты недостоверные, то и искомые величины будут недостоверными.

Слежу за темой, но вижу ты уходишь в сторону. Проверь еще раз свои матрицы, правильно их приведи по Гаусу, найди обратные итд...
Не возможно привести матрицу, и оставить диагональ не из единиц.
В этом сам алгоритм.

А по поводу метода наименьших квадратов - какие проблемы?
http://ru.wikipedia.org/wiki/Метод_наименьших_квадратов

[Сергей Дубина]

Сообщение от Елпанов Евгений А по поводу метода наименьших квадратов - какие проблемы?

На странице рассмотрен пример. Проверил ход решения результат верный, но подставив полученные значения получим:
16 --> 18,62177758
32 --> 28,51259293
55 --> 55,9385323
32 --> 31,19376103
29 --> 29,7811114
для:
x = 3,55
y = − 0,109
Так и должно быть?

[Елпанов Евгений]

Сообщение от Сергей Дубина Так и должно быть?

Да! Так и должно быть!

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

[Сергей Дубина]

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

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

[Елпанов Евгений]

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

[Сергей Дубина]

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

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

[Елпанов Евгений]

посмотри здесь:
наглядное пояснение смысла метода на примере точек и линии...

Для личной беседы - пиши в аську.

[Сергей Дубина]

Сообщение от Елпанов Евгений посмотри здесь:

С апроксимацией и МНК более менее понятно, а вот насчет моей проблемы (см. пост #11), не могли бы направить в нужное русло.

[Елпанов Евгений]

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

так ведь любое, из бесконечного количества решений, решение - то, что вам надо...
Зачем здесь метод наименьших квадратов?
Может просто методом Гауса решить систему? Лишние уравнения просто отбросив.
Смысл метода наименьших квадратов - решить систему, у которой нет единого решения. На примере точек и одной прямой - любое из уравнений матрицы - прямая, но все уравнения из матрицы описывают разные прямые. В этом случае, если требуется найти хоть какое решение, которое будет в равной степени близко ко всем описанным прямым и используют метод наименьших квадратов.
т.е. если есть облако точек и необходимо его описать одной точкой -
1. можно взят среднюю точку по массе - все точки имеют одинаковый вес.
2. можно взять середину описанной окружности - как ни крути, эта середина будет в одном и том же месте.
3. можно взять за основу, произвольную точку из этого облака - точка всегда внутри облака.

То же самое и с методом наименьших квадратов. Допустим есть список точек после сканирования плоскости и необходимо найти эту плоскость.
Допустим точность сканера не высока и мы имеем десяток точек на стене и нам необходимо ее отрисовать плоскостью, но в действительности, нет ни одной четверки точек лежащих на одной плоскости. Как быть? Именно для этого случая и применяют метод. Можно найти среднее решение, которое не обязательно совпадет хотя бы с одним из уравнений, но будет описывать некоторое среднее значение.

(Нужен алгоритм, умножаем, транспонируем .... получаем произвольное решение)

Смотри сообщение №6. Там ссылка на простое описание алгоритма метода наименьших квадратов, решение самой матрицы - метод Гаусса

[Сергей Дубина]

Вопрос с МНК, апроксимацией и методом Гаусса для меня исчерпан, и заинтересовался им "на песпективу".
Я решаю другую задачу. Мне просто не хочется создавать новую тему (эта "пылится" уже полгода и каким то чудом ожила).
Если есть желание могу вынести мой субвопрос в отдельную тему.

[Елпанов Евгений]

Код:

[Выделить все]

(defun c:cce (/ LP LST)
 ;; by ElpanovEvgeniy
 ;; Calculation of the connected equations
 ;; 13.11.2006
 ;; (c:cce)
 (if (setq lst (text-equation-list))
  (progn
   (setq lp  (append
              (ACAD_STRLSORT
               (rem-dupl (vl-remove "=" (apply (function append) (mapcar (function car) lst))))
              ) ;_  ACAD_STRLSORT
              '("=")
             ) ;_  append
         lst (vl-remove-if
              (function (lambda (x) (vl-every (function zerop) x)))
              (gauss (mapcar (function (lambda (x)
                                        (mapcar (function (lambda (a)
                                                           (if (assoc a (cadr x))
                                                            (cdr (assoc a (cadr x)))
                                                            0.
                                                           ) ;_  if
                                                          ) ;_  lambda
                                                ) ;_  function
                                                lp
                                        ) ;_  list
                                       ) ;_  lambda
                             ) ;_  function
                             lst
                     ) ;_  mapcar
              ) ;_  gauss
             ) ;_  vl-remove-if
   ) ;_  setq
   (if (= (1- (length (car lst))) (length lst))
    (if (and (= (length (vl-remove-if (function zerop) (last lst))) 1)
             (not (zerop (last (last lst))))
        ) ;_  and
     (entmakex (list '(0 . "MTEXT")
                     '(100 . "AcDbEntity")
                     '(100 . "AcDbMText")
                     (cons 10 (getpoint "\nselect a point to place the result:"))
                     '(1 . "This set of equations\n has no solutions!")
               ) ;_  list
     ) ;_  entmakex
     (entmakex
      (list
       '(0 . "MTEXT")
       '(100 . "AcDbEntity")
       '(100 . "AcDbMText")
       (cons 10 (getpoint "\nselect a point to place the result:"))
       (cons
        1
        (apply
         (function strcat)
         (mapcar
          (function
           (lambda (a b)
            (strcat
             a
             "="
             (rtos (apply (function /) (reverse (vl-remove-if (function zerop) b))) 2 8)
             "\n"
            ) ;_  strcat
           ) ;_  lambda
          ) ;_  function
          lp
          (mapcar
           (function (lambda (x) (cdr (reverse (cons (last x) x)))))
           (reverse
            (gauss (reverse
                    (mapcar (function (lambda (x) (cdr (reverse (cons (last x) x))))) lst)
                   ) ;_  reverse
            ) ;_  gauss
           ) ;_  reverse
          ) ;_  reverse
         ) ;_  mapcar
        ) ;_  apply
       ) ;_  cons
      ) ;_  list
     ) ;_  entmakex
    ) ;_  if
    (progn
     (setq
      lst (reverse
           (mapcar
            (function (lambda (x)
                       (vl-remove-if (function (lambda (x1) (zerop (cdr x1))))
                                     (mapcar (function cons) lp x)
                       ) ;_  vl-remove-if
                      ) ;_  lambda
            ) ;_  function
            (mapcar
             (function (lambda (x) (cdr (reverse (cons (last x) x)))))
             (reverse
              (gauss
               (reverse
                (mapcar (function (lambda (x) (cdr (reverse (cons (last x) x))))) lst) ;_  mapcar
               ) ;_  reverse
              ) ;_  gauss
             ) ;_  reverse
            ) ;_  mapcar
           ) ;_  mapcar
          ) ;_  reverse
     ) ;_  setq
     (if (= 1 (length (last lst)))
      (entmakex (list '(0 . "MTEXT")
                      '(100 . "AcDbEntity")
                      '(100 . "AcDbMText")
                      (cons 10 (getpoint "\nselect a point to place the result:"))
                      '(1 . "This set of equations\n has no solutions!")
                ) ;_  list
      ) ;_  entmakex
      (entmakex
       (list
        '(0 . "MTEXT")
        '(100 . "AcDbEntity")
        '(100 . "AcDbMText")
        (cons 10 (getpoint "\nselect a point to place the result:"))
        (cons
         1
         (apply
          (function strcat)
          (mapcar
           (function
            (lambda (x)
             (apply
              (function strcat)
              (append
               (mapcar
                (function (lambda (a)
                           (if (= (car a) "=")
                            (strcat (car a) (rtos (cdr a) 2 8))
                            (cond ((equal (cdr a) 1. 1e-8) (car a))
                                  ((equal (cdr a) -1. 1e-8) (strcat "-" (car a)))
                                  ((minusp (cdr a)) (strcat (rtos (cdr a) 2 8) (car a)))
                                  (t (strcat "+" (rtos (cdr a) 2 8) (car a)))
                            ) ;_  cond
                           ) ;_  if
                          ) ;_  lambda
                ) ;_  function
                x
               ) ;_  mapcar
               '("\n")
              ) ;_  append
             ) ;_  apply
            ) ;_  lambda
           ) ;_  function
           (reverse
            (Simplification-equations
             (mapcar
              (function (lambda (x)
                         (append (list (car x) (last x))
                                 (mapcar (function (lambda (a) (cons (car a) (- (cdr a)))))
                                         (reverse (cdr (reverse (cdr x))))
                                 ) ;_  mapcar
                         ) ;_  append
                        ) ;_  lambda
              ) ;_  function
              lst
             ) ;_  mapcar
            ) ;_  Simplification-equations
           ) ;_  reverse
          ) ;_  mapcar
         ) ;_  apply
        ) ;_  cons
       ) ;_  list
      ) ;_  entmakex
     ) ;_  if
    ) ;_  progn
   ) ;_  if
  ) ;_  progn
 ) ;_  if
) ;_  defun
(defun Simplification-equations (l)
 ;; By ElpanovEvgeniy
 ;; Simplification of the connected equations
 ;|
(Simplification-equations '((("C" . 1.0) ("=" . -3.0) ("D" . -1.0))
       (("B" . 1.0) ("=" . -3.0) ("C" . 2.0))
       (("A" .  1.0) ("=" . 4.0) ("B" .  2.0))
      )
)
  ;; =>
'((("C" . 1.0) ("=" . -3.0) ("D" . -1.0))
  (("B" . 1.0) ("=" . -9.0) ("D" . -2.0))
  (("A" . 1.0)  ("=" . -14.0) ("D" .  -4.0))
 )
  |;
 (if (cdr l)
  (cons
   (cons (cons (caaar l) 1.)
         (mapcar (function (lambda (x) (cons (car x) (/ (cdr x) (cdaar l))))) (cdar l))
   ) ;_  cons
   (Simplification-equations
    (cons
     (cons
      (cons (caaadr l) 1.)
      (cons
       (cons
        "="
        (+ (cdr (assoc "=" (cadr l)))
           (* (/ (cdr (assoc (caaar l) (cadr l))) (cdaar l)) (cdr (assoc "=" (car l))))
        ) ;_  +
       ) ;_  cons
       (mapcar (function
                (lambda (x)
                 (cons (car x) (* (cdr x) (/ (cdr (assoc (caaar l) (cadr l))) (cdaar l))))
                ) ;_  lambda
               ) ;_  function
               (cddar l)
       ) ;_  mapcar
      ) ;_  cons
     ) ;_  cons
     (cddr l)
    ) ;_  cons
   ) ;_  Simplification-equations
  ) ;_  cons
  (list
   (cons (cons (caaar l) 1.)
         (mapcar (function (lambda (x) (cons (car x) (/ (cdr x) (cdaar l))))) (cdar l))
   ) ;_  cons
  ) ;_  list
 ) ;_  if
) ;_  defun
(defun rem-dupl (lst)
 ;; by ElpanovEvgeniy
 ;; Deleting of all identical objects
 ;; (rem-dupl '(1 2 3 2 1))=>(1 2 3)
 (if lst
  (cons (car lst) (rem-dupl (vl-remove (car lst) (cdr lst))))
 ) ;_  if
) ;_  defun
(defun gauss (lst)
 ;; by ElpanovEvgeniy
 ;; Implementation Gaussian elimination
 ;; For text:
                ;  1x+2y+3z=2
                ;  10x+1y+8z=17
                ;  7z+2y=5
 ;; (gauss '((1.0 2.0 3.0 2.0) (10.0 1.0 8.0 17.0) (0.0 2.0 7.0 5.0)))
 ;; =>
 ;; ((1.0 2.0 3.0 2.0) (0.0 -19.0 -22.0 -3.0) (0.0 0.0 4.68421 4.68421))
 ;;(gauss lst)
 (if (car lst)
  (if (zerop (caar lst))
   (if (vl-every (function zerop) (mapcar (function car) lst))
    (if (cdr lst)
     (cons (car lst)
           (mapcar (function (lambda (x) (cons 0. x)))
                   (gauss (mapcar (function cdr) (cdr lst)))
           ) ;_  mapcar
     ) ;_  cons
     lst
    ) ;_  if
    (gauss
     (cons (mapcar (function +)
                   (car lst)
                   (car (vl-remove-if (function (lambda (x) (zerop (car x)))) (cdr lst)))
           ) ;_  mapcar
           (cdr lst)
     ) ;_  cons
    ) ;_  gauss
   ) ;_  if
   (cons (car lst)
         (mapcar
          (function (lambda (x) (cons 0. x)))
          (gauss
           (mapcar (function (lambda (x / i)
                              (setq i (/ (car x) (caar lst)))
                              (mapcar (function -)
                                      (cdr x)
                                      (mapcar (function (lambda (x1) (* x1 i))) (cdar lst))
                              ) ;_  mapcar
                             ) ;_  lambda
                   ) ;_  function
                   (cdr lst)
           ) ;_  mapcar
          ) ;_  test
         ) ;_  mapcar
   ) ;_  cons
  ) ;_  if
 ) ;_  if
) ;_  defun

(defun text-equation-list (/ s)
 ;; By ElpanovEvgeniy
 ;; select text equation and create list
 ;; For text:
                ;  1x+2y+3z=2
                ;  10x+1y+8z=17
                ;  7z+2y=5
 ;; (text-equation-list)
 ;; =>
 ;|'((("X" "Y" "Z" "=") (("X" . 1.0) ("Y" . 2.0) ("Z" . 3.0) ("=" . 2.0)))
       (("X" "Y" "Z" "=") (("X" . 10.0) ("Y" . 1.0) ("Z" . 8.0) ("=" . 17.0)))
       (("Z" "Y" "=") (("Z" . 7.0) ("Y" . 2.0) ("=" . 5.0)))
     )|;
 (if (setq s (ssget '((0 . "TEXT"))))
  (mapcar
   (function (lambda (str)
              (setq str
                    (mapcar
                     (function cons)
                     (mapcar
                      (function VL-PRINC-TO-STRING)
                      (read (strcat "("
                                    (VL-STRING-TRANSLATE "0123456789.=+-" "              " str) ;_  VL-STRING-TRANSLATE
                                    "=)"
                            ) ;_  strcat
                      ) ;_  read
                     ) ;_  mapcar
                     (mapcar
                      (function FLOAT)
                      (subst
                       -1
                       '-
                       (subst
                        1
                        '+
                        (read
                         (strcat "("
                                 (VL-STRING-TRANSLATE
                                  "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ="
                                  "                                                     "
                                  (cond ((WCMATCH (substr str 1 1) "[a-z],[A-Z]") (strcat "1" str))
                                        ((WCMATCH (substr str 1 2) "-[a-z],-[A-Z]")
                                         (strcat "-1" (substr str 2))
                                        )
                                        (t str)
                                  ) ;_  cond
                                 ) ;_  VL-STRING-TRANSLATE
                                 ")"
                         ) ;_  strcat
                        ) ;_  read
                       ) ;_  subst
                      ) ;_  subst
                     ) ;_  mapcar
                    ) ;_  mapcar
              ) ;_  setq
              (list (mapcar (function car) str) str)
             ) ;_  lambda
   ) ;_  function
   (mapcar (function (lambda (x) (cdr (assoc 1 (entget x)))))
           (vl-sort
            (vl-remove-if (function listp) (mapcar (function cadr) (ssnamex s)))
            (function
             (lambda (a b) (> (caddr (assoc 10 (entget a))) (caddr (assoc 10 (entget b)))))
            ) ;_  function
           ) ;_  vl-sort
   ) ;_  mapcar
  ) ;_  mapcar
 ) ;_  if
) ;_  defun

[Сергей Дубина]

Попробую переформатировать вопрос:
Имееется прямоугольная матрица А (n*m; где n<m) и матрица-столбец B (1*m).
Необходимо найти матрицу-столбец Х (1*n) (любое правильное решение из бесконечного количества)
n - количество неизвестных
m - количество уравнений
нужно решение в матричном виде (например Х={A}*{A}т........)

[slogos]

Решение такое. 1)По исходной матрице А нужно составить матрицу нормальных уравнений B, в которой коэффициентами будут производные по соответствующим переменным (т.н. матрица Якоби) 2)а далее все просто X = B(обр)*L. где L - столбец (правая часть) системы. Он получается одновременно с нормальной матрицей (как расширенной)
Вот и вся теория.

[Сергей Дубина]

X=(((АтА)-1)Ат)B

[slogos]

Совершенно верно, но уже (АтА) с большими коэффициентами, ((АтА)-1) - с очень большими, Е = (АтА)*((АтА)-1)) - "не такаю уж еденичная" и в конечном итоге X = (x1, x2, ...Xn) - ерунда! Вот в чем фикус.