Matlab. Функция для обращения плохо обусловленной матрицы

[Simply me]

Здравствуйте! Помогите, пожалуйста, найти или написать функцию для обращения плохо обусловленных матриц в MATLAB. Из встроенных матлабовских пробовала функции inv и pinv (псевдообратная с использованием svd-разложения, 32-битная). Проверяю на матрице [1 2 3; 2 4 6; 1 1 1]. Не работают. Реализовала через нахождение LU-разложения. Ведущий элемент пробовала выбирать как максимальный и как элемент из строки, которая образует матрицу с наименьшим числом обусловленности. Для плохо обусловленных не работает. Ещё на файлообменнике Матлаб нашла функции inverse, inv_chol, ginv, qrginv, pseudoinverse. Тоже плохая точность. Как можно улучшить точность?

[one man]

DIGITS

Пользуетесь?

[мат-ламер]

Simply me в сообщении #1011744 писал(а):

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

Simply me в сообщении #1011744 писал(а):

Проверяю на матрице [1 2 3; 2 4 6; 1 1 1].

А что, это плохо обусловленная матрица?

-- Ср май 06, 2015 22:26:46 --

one man в сообщении #1011797 писал(а):

DIGITS

Пользуетесь?

И какую точность надо задать для обращения такой матрицы?

-- Ср май 06, 2015 22:27:29 --

А псевдообратная для такой матрицы должна получиться.

[one man]

Был околдован словами “плохо обусловленных” и даже не глянул на пример.

[Simply me]

Похоже, путаю плохо обусловленные матрицы и матрицы с линейно зависимыми строками (т.е. когда ранг меньше числа строк). Хорошо бы, чтобы и для тех, и для других работало обращение. Для A=[1 2 3; 2 4 6; 1 1 1] А*pinv(A)*A-A получаются числа порядка 10е-15. Нормально, наверное.
Вообще у меня в программе есть обращение плохо обусловленной матрицы P. Использовала pinv. На 64 разрядном комп-е было нормально. На 32-разрядном запустила, получается плохо. Тогда занулила начальные условия (матрица P зависит от начальных условий), запустила на 64, получиласьплохая точность. Слышала, что у pinv 32-битная точность. Может, у нее есть какой-нибудь 64-битный аналог?

[Geen]

Simply me в сообщении #1012038 писал(а):

Слышала, что у pinv 32-битная точность.

Это как?? Matlab не работает с числами одинарной точности по умолчанию, даже на 16-битной архитектуре.

[Simply me]

Не знаю, может, перепутала. Все равно pinv не работает для плохо обусловленных матриц. Пробовала передать второй параметр: pinv(P, tol). Если так вызывать, то все сингулярные числа, которые меньше tol, будут зануляться. Задала tol меньше, чем был по умолчанию. Но все равно не помогло. Ещё пробовала метод регуляризации Тихонова. Получается не лучше, чем pinv (а если взять параметр больше 1, то намного хуже). Еще пробовала hpf (high precision floating point). Зависает. Может, можно ещё что-нибудь придумать?

[Geen]

Simply me в сообщении #1012038 писал(а):

получаются числа порядка 10е-15.

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

[Simply me]

Такая точность 10е-15 получается для матрицы A=[1 2 3; 2 4 6; 1 1 1]. Я её сама придумала А для той матрицы, которая в программе, плохо получается. На 64 с ненулевыми условиями хорошо, на 64 с нулевыми плохо, на 32 с ненулевыми плохо. Как сделать, чтобы и для плохо обусловленных хорошо получалось?

[Simply me]

Помогите, пожалуйста! Наверняка кто-нибудь сталкивался с такой проблемой, потому что когда имеешь дело с реальными данными не бывает все гладко, например, попадаются плохо обусловленные матрицы. Может, есть в матлаб функция для обращения, более точная, чем pinv? Или, может, есть какой-нибудь алгоритм, устойчивый к ошибкам округления? Сама пробовала закодить метод регуляризации и нахождение обратной через lu-факторизацию. Ведущий элемент в lu-факторизации пробовала выбирать несколькими способами. Ещё нашла одну статью. Там как раз про недостатки функции pinv. Предлагается использовать QR-разложение и приводится код функции qrginv. Для моей матрицы это все не работает.

[Alexu007]

А что вы (или программа) проделываете с целыми числами A=[1 2 3; 2 4 6; 1 1 1], что получается потеря точности?

[Geen]

Simply me в сообщении #1012358 писал(а):

Для моей матрицы это все не работает.

Загадочная матрица.
А зачем Вам вообще псевдообратная, и именно от этой матрицы? Может попробовать переформулировать задачу? Иногда помогает.

[Simply me]

Geen, задача идентификации. То есть, есть входной сигнал u и выходной сигнал y. Нужно восстановить параметры системы, потом подставить в эту систему найденные параметры и входной сигнал, и полученный выходной сигнал y' должен быть похожим на y. Система диф. уравнений. Она преобразуется к системе линейных уравнений: y=P*x. Нужно найти x, а для этого нужно найти обратную к P. Раньше все везде работало, графики y и y' бились. А потом решила сделать, чтобы алгоритм был устойчивым к выбросам. Алгоритм немного переделала, pinv стала вызываться в цикле. Входные данные изменила: занулила 20% массива y в середине. Работало. Запустила на 32 битном компьютере-не работает. Задала у системы нулевые начальные условия на 64 битном - не работает. Вообщем опытным путём установлено, что когда матрица P плохо обусловлена, тогда и не работает. Нужно, чтобы было точнее. Матрица P у меня есть в матлабовском файле *.mat, но сюда по-моему нельзя прикреплять файлы.

[Geen]

Simply me в сообщении #1012387 писал(а):

А потом решила сделать, чтобы алгоритм был устойчивым к выбросам.

В этом ключевая ошибка - по своему опыту сужу

Я правильно понимаю, что у $P$ нулевой детерминант? (она вообще квадратная?) Она фиксированная?
Что такое выбросы?
Почему бы не попробовать просто минимизировать невязку?

[Simply me]

P не квадратная. В некоторые функции, в частности в pinv, передаю P, а в те, которые принимают только квадратные матрицы, передают P'*P. У P'*P нулевой детерминант. У P большое число обусловленности 10е+21. У P'*P ещё больше 10е+35 (для сравнения если с ненулевыми условиями cond(P)=10e+8 и pinv(P) нормально работает). P не совсем фиксированная, потому что она получается через u и y, y получается через u, а u через функцию randn. Но в принципе P можно считать постоянной при каждом запуске, не сильно меняется. Я вообще её сохранила в файл и с ней экспериментирую.
Выбросы - это отклонения в исходных данных. Их нужно учитывать. Потому что модель будет пытаться как бы "угодить" всем замерам, в том числе выбросам, и в угоду им сдвинется. Я из уравнения y=P*x нахожу x методом наим квадратов. Потом нахожу невязку, сортирую, нахожу индексы 20% уравнений, у которых наиб невязка, отбрасываю эти индексы из y и P. Снова решаю y=P*x. Полученный x подставляю уже в первоначальную систему y=P*x. Опять считаю невязку и т.д. В результате для работающего случая (64 бит, ненулевые условия) графики y и y' полностью совпадают ( на самом деле не совсем, но визуально кажется, что полностью) везде, кроме 20% замеров в серединке. В неработающих случаях они в зоне выбросов не совпадают, но и до и после выбросов тоже не очень совпадают. Если улучшить точность вычислений, то это все исправится.