2014 dxdy logo

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

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




Начать новую тему Ответить на тему На страницу 1, 2  След.
 
 Обращение и псевдо обращение матриц
Сообщение04.07.2012, 12:03 
Аватара пользователя


31/10/08
1244
Какие есть алгоритмы для псевдообращений матриц? Существует ли блочный алгоритм?
Я знаю только один алгоритм - это SVD.

 Профиль  
                  
 
 Re: Обращение и псевдо обращение матриц
Сообщение04.07.2012, 13:48 
Заслуженный участник


11/05/08
31922
SVD -- штука затратная в том смысле, что не реализуется за конечное количество шагов. Псевдообращение же само по себе -- вполне реализуемо.

Псевдообратная матрица ставит в соответствие каждому вектору $f$ правых частей нормальное псевдорешение, т.е. решение $x$ системы $A^*Ax=A^*f$, имеющее минимальную норму. Собственно решение можно найти, скажем, обычным методом Гаусса. Причём если решений много, то достаточно найти одно; тогда нормальное решение получается просто ортогонализацией этого одного решения ядру матрицы $A^*A$.

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

 Профиль  
                  
 
 Re: Обращение и псевдо обращение матриц
Сообщение04.07.2012, 15:26 
Аватара пользователя


31/10/08
1244
ewert
Цитата:
SVD -- штука затратная в том смысле, что не реализуется за конечное количество шагов.

А то в книге говорится Уилкинсон_Дж.,_Райнш_Ц.-Справочник_алгоритмов_на_языке_АЛГОЛ._Линейная_алгебра(1976)
что по расчётам в среднем 2 итерации. И в общем случае там предложено ограничить 32 интерациями на сингулярное число.
Во вторых как я понимаю проблемы сходимости можно решить и свести к числу разрядов используемых для представления чисел. Процесс с математической точкой может быть и бесконечным, но алгоритмически он ограничен машинной точностью(машинной ошибкой).

Цитата:
Псевдообратная матрица ставит в соответствие каждому вектору $f$ правых частей нормальное псевдорешение, т.е. решение $x$ системы $A^*Ax=A^*f$, имеющее минимальную норму.
Совсем не понял фразы. :-( Что такое $f$ и $x$? что такое нормально решение?

Цитата:
какие подходы тут считаются наиболее модными.

Мне модный не нужен. Мне нужен универсальный метод, а именно чтобы работал с действительной матрице, а лучше комплексной(на практике лучше 2 метода, а то комплексные числа много памяти едят).
Чтобы работал не только с квадратной матрицей, но и с прямоугольной. Что-бы не было проблем со сходимостью. Что-бы не было проблем с "обусловленностью". Не было проблем с ошибками округления.

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

У меня на практике порог не спасает и я не понимаю почему.

 Профиль  
                  
 
 Re: Обращение и псевдо обращение матриц
Сообщение04.07.2012, 15:59 
Заслуженный участник


11/05/08
31922
Pavia в сообщении #592048 писал(а):
Что такое $f$ и $x$? что такое нормально решение?

Псевдорешением системы $Ax=f$ называется (в конце концов) решение системы $A^*Ax=A^*f$.

Нормальным решением системы $Bx=f$ называется решение, имеющее наименьшую норму (если, конечно, решений бесконечно много). Или, что эквивалентно -- решение, ортогональное ядру матрицы.

Нормальным псевдорешением системы $Ax=f$ называется её псевдорешение, являющееся нормальным.

Псевдообратной к матрице $A$ называется матрица $A^+$, которая каждый вектор $f$ переводит в вектор $x=A^+f$, являющийся нормальным псевдорешением системы $Ax=f$.

Все эти вещи имеют смысл для матриц любого размера (т.е. не обязательно квадратных) и легко реализуются за конечное количество шагов, так что вопрос о сходимости просто не встаёт. Есть ещё явная система уравнений Пенроуза для псевдообратной матрицы $A^+$, но она какая-то занудная.

Pavia в сообщении #592048 писал(а):
Что-бы не было проблем с "обусловленностью". Не было проблем с ошибками округления.

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

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

 Профиль  
                  
 
 Re: Обращение и псевдо обращение матриц
Сообщение06.07.2012, 09:24 
Аватара пользователя


06/07/12
70
Pavia в сообщении #591971 писал(а):
Какие есть алгоритмы для псевдообращений матриц? Существует ли блочный алгоритм?
Я знаю только один алгоритм - это SVD.

SVD - единственный алгоритм для общего случая (Мура-Пенроуза матрица).
Блочные алгоритмы существуют, например, в пакете LAPACK. Его достоинством является то, что он во-первых создан при содействии известных американских математиков и во-вторых постоянно обновляется. В-третьих, его одинаково легко использовать как с двойной, так и с четверной (IEEE 754-2008: скорость расчета на порядок падает) точностью. Но для больших матриц лучше использовать профессиональную реализацию пакета LAPACK - Intel MKL: скорость счета увеличивается на порядок (правда, четверной точности там нет, но под AVX скорость возрастает почти в 20 раз). Под Linux некоммерческая версия бесплатна.

 Профиль  
                  
 
 Re: Обращение и псевдо обращение матриц
Сообщение07.07.2012, 16:05 
Заслуженный участник


11/05/08
31922
Ну пожалуйста. Сначала программируем нахождение псевдорешения (какого-либо) системы $Ax=f$:

код: [ скачать ] [ спрятать ]
Используется синтаксис Matlab M
function [z,x]=psol(A,f)
% Псевдорешение системы Ax=f
% (метод Гаусса с полным выбором главного элемента)
% z - ядро матрицы (или пусто)

tol=1e-10;
x=[];
z=[];

[m,n]=size(A);
[mr,nr]=size(f);
if mr~=m
        disp('Несогласованность размеров')
        return
end
       
B=A'*A;
mb=max(max(abs(B)));    %  Оценка нормы матрицы
B=[[B, A'*f]; 1:(n+nr)];    %  Расширенная матрица с добавленной строкой перестановок

nz=0;    %  Для размерности ядра
for i=1:n
    [mc,ir]=max(abs(B(i:n,i:n)));
    [mm,ic]=max(mc);
    ir=ir(ic)+i-1;
    ic=ic+i-1;      %  (ir,ic) - индексы главного элемента
    if mm<mb*tol
        nz=n-i+1;
        break;     %  Метод Гаусса завершен - остальные строки нулевые
    end
    if mb<mm,    mb=mm;    end
    rr=B(i,:);    B(i,:)=B(ir,:);    B(ir,:)=rr;    %  Переставили строки
    rr=B(:,i);    B(:,i)=B(:,ic);    B(:,ic)=rr;    %  Переставили столбцы
    B(i,:)=B(i,:) / B(i,i);
    B([1:i-1,i+1:n], :)=B([1:i-1,i+1:n], :) - B([1:i-1,i+1:n], i) * B(i,:);    %  Очистили i-й столбец
end

ix=B(n+1, 1:n-nz);    %  Индексы связанных переменных
iz=B(n+1, n-nz+1:n);    %  Индексы свободных переменных

x=zeros(n, nr);
x(ix, :)=B(1:n-nz, n+1:n+nr);    %  Частное решение

if nz>0        %  Базис ядра:
    z=eye(n,nz);
    z(iz, :)=z(1:nz, :);
    z(ix, 1:nz)=-B(1:n-nz, n-nz+1:n);
end

На выходе x -- это матрица, каждый столбец которой есть псевдорешение для соответствующего столбца матрицы f; матрица z (если не пустая) состоит из столбцов, образующих базис ядра исходной матрицы. Теперь, чтобы найти псевдообратную матрицу, достаточно запустить эту процедуру с единичной матрицей в качестве f, после чего ортогонализовать каждый столбец полученной x к ядру:

код: [ скачать ] [ спрятать ]
Используется синтаксис Matlab M
function x=pseudo(A)
% Нахождение псевдообратной матрицы
% (на основе метода Гаусса с полным выбором главного элемента)

[m,n]=size(A);
f=eye(m,m);

[z,x]=psol(A,f);
[i,nx]=size(x);
[i,nz]=size(z);

for i=1:nz        %  Ортонормируем ядро по Граму-Шмидту
    for k=1:i-1
        z(:, i) = z(:, i) - z(:, k) * sum(z(:, i).*z(:, k));
    end
    z(:, i) = z(:, i) / sqrt(sum(z(:, i).^2));
end

for i=1:nx        %  Ортогонализуем псевдорешения к ядру
    for k=1:nz
            x(:, i) = x(:, i) - z(:, k) * sum(x(:, i).*z(:, k));
    end
end

Вот тестовая программка:

код: [ скачать ] [ спрятать ]
Используется синтаксис Matlab M
m=6;        %  Высота матрицы
n=3;        %  Ширина матрицы
d=1;        %  Неполнота ранга матрицы

aa=rand(m, n);
a=aa*0;

if m>n
    c=rand(n, n-d);
    for i=1:n
        for k=1:n-d
            a(:,i)=a(:,i) + aa(:,k) * c(i,k);
        end
    end
else
    c=rand(m, m-d);
    for i=1:m
        for k=1:m-d
            a(i,:)=a(i,:) + aa(k,:) * c(i,k);
        end
    end
end

a

tic, b=pseudo(a), toc
tic, c=pinv(a), toc

max(max(abs(b-c)))

%max(max(abs(b-a^(-1))))
%max(max(abs(c-a^(-1))))

В принципе, при любом варианте исходной матрицы независимо от её размеров и неполноты ранга получается то же, что и при использовании стандартной матлабовской функции pinv. Правда, точность на пару знаков ниже, но я ведь и не пытался выжать все соки. Тем более не пытался что-то оптимизировать.

 Профиль  
                  
 
 Re: Обращение и псевдо обращение матриц
Сообщение07.07.2012, 17:26 
Аватара пользователя


31/10/08
1244
ewert
Спасибо за код. Буду разбираться.

Матрица у меня плохо обусловлена.
Код:
a=[4097 -3372.3292364 5730820823.8 -14151482820 1.4429140017e+16 -5.9384620888e+16 4.3249891776e+22; -3372.3292364 5730820823.8 -14151482820 1.4429140017e+16 -5.9384620888e+16 4.3249891776e+22 -2.4919888057e+23; 5730820823.8 -14151482820 1.4429140017e+16 -5.9384620888e+16 4.3249891776e+22 -2.4919888057e+23 1.4116057367e+29; -14151482820 1.4429140017e+16 -5.9384620888e+16 4.3249891776e+22 -2.4919888057e+23 1.4116057367e+29 -1.0457268159e+30; 1.4429140017e+16 -5.9384620888e+16 4.3249891776e+22 -2.4919888057e+23 1.4116057367e+29 -1.0457268159e+30 4.8465906922e+35; -5.9384620888e+16 4.3249891776e+22 -2.4919888057e+23 1.4116057367e+29 -1.0457268159e+30 4.8465906922e+35 -4.3882409082e+36; 4.3249891776e+22 -2.4919888057e+23 1.4116057367e+29 -1.0457268159e+30 4.8465906922e+35 -4.3882409082e+36 1.7209134454e+42]

pseudo(a) работает ещё хуже чем pinv(a).
Матрица квадратная и ту можно использовать inv. inv(a) лучше работает чем pinv(a)

Нашел такую работу http://www.ict.edu.ru/vconf/files/11475.pdf
Думаю по пробовать алгоритм Эрмита, шестое чувство говорит что он должен лучше работать.

А вот точное решение с использованием рациональных чисел заинтриговала. Стоит ли прорабатывать или дальше с плавающими цифрами воевать?

 Профиль  
                  
 
 Re: Обращение и псевдо обращение матриц
Сообщение07.07.2012, 17:55 
Аватара пользователя


06/07/12
70
Pavia,

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

 Профиль  
                  
 
 Re: Обращение и псевдо обращение матриц
Сообщение07.07.2012, 18:29 
Аватара пользователя


31/10/08
1244
ogaman
В матлабе pinv реализован через svd.
У pinv(a) проблема с погрешностью растёт экспоненциально, с ростом матрицы(в моих задачах). Что не в какие ворота не лезет.
Цитата:
реализованный профессионалами мирового уровня?

Не боги горшки обжигают.

 Профиль  
                  
 
 Re: Обращение и псевдо обращение матриц
Сообщение07.07.2012, 18:53 
Аватара пользователя


06/07/12
70
Pavia,

лучше SVD пока ничего не придумали. И Вам надо использовать не pinv(a), а pinv(a,tol). Кроме того я отмечал, что решение с помощью рекомендованного мною кода можно находить с четверной точностью: в Вашем случае это существенно.

 Профиль  
                  
 
 Re: Обращение и псевдо обращение матриц
Сообщение08.07.2012, 13:58 
Аватара пользователя


06/07/12
70
Хотелось бы отметить, что новейшие реализации SVD довольно сложны: например, самый быстрый алгоритм, реализованный на языке высокого уровня, содержит около 80 процедур и функций и занимет на диске около одного мегабайта. А на реализацию некоторых частей этого алгоритма ушло не менее 10-ти лет: т.е алгоритм был разработан, а его программная реализация потребовала 10-летней работы.

 Профиль  
                  
 
 Re: Обращение и псевдо обращение матриц
Сообщение08.07.2012, 14:54 
Заслуженный участник


11/05/08
31922
Pavia в сообщении #593104 писал(а):
pseudo(a) работает ещё хуже чем pinv(a).

На самом деле в определённом смысле наоборот. Вот проверка уравнений Пенроуза для матрицы, полученной двумя способами:

код: [ скачать ] [ спрятать ]
Используется синтаксис Matlab M
a=[ 4097             -3372.3292364      5730820823.8     -14151482820       1.4429140017e+16 -5.9384620888e+16  4.3249891776e+22
   -3372.3292364      5730820823.8     -14151482820       1.4429140017e+16 -5.9384620888e+16  4.3249891776e+22 -2.4919888057e+23
    5730820823.8     -14151482820       1.4429140017e+16 -5.9384620888e+16  4.3249891776e+22 -2.4919888057e+23  1.4116057367e+29      
   -14151482820       1.4429140017e+16 -5.9384620888e+16  4.3249891776e+22 -2.4919888057e+23  1.4116057367e+29 -1.0457268159e+30
    1.4429140017e+16 -5.9384620888e+16  4.3249891776e+22 -2.4919888057e+23  1.4116057367e+29 -1.0457268159e+30  4.8465906922e+35
   -5.9384620888e+16  4.3249891776e+22 -2.4919888057e+23  1.4116057367e+29 -1.0457268159e+30  4.8465906922e+35 -4.3882409082e+36
    4.3249891776e+22 -2.4919888057e+23  1.4116057367e+29 -1.0457268159e+30  4.8465906922e+35 -4.3882409082e+36  1.7209134454e+42];

a=a*1e-42;

disp('B=pseudo(A)')
b=pseudo(a);
disp('|ABA-A|=')
max(max(abs(a*b*a-a)))
max(max(abs((a*b-eye(7))*a)))
max(max(abs(a*(b*a-eye(7)))))
disp('|BAB-B|=')
max(max(abs(b*a*b-b)))
max(max(abs(b*(a*b-eye(7)))))
max(max(abs((b*a-eye(7))*b)))
disp('|A''B''-BA|=')
max(max(abs(a'*b'-b*a)))
disp('|B''A''-AB|=')
max(max(abs(b'*a'-a*b)))

disp('================================')

disp('C=pinv(A)')
c=pinv(a);
disp('|ACA-A|=')
max(max(abs(a*c*a-a)))
max(max(abs((a*c-eye(7))*a)))
max(max(abs(a*(c*a-eye(7)))))
disp('|CAC-C|=')
max(max(abs(c*a*c-c)))
max(max(abs(c*(a*c-eye(7)))))
max(max(abs((c*a-eye(7))*c)))
disp('|A''C''-CA|=')
max(max(abs(a'*c'-c*a)))
disp('|C''A''-AC|=')
max(max(abs(c'*a'-a*c)))

disp('================================')

disp('|B|=')
max(max(abs(b)))
disp('|C|=')
max(max(abs(c)))
disp('|B-C|=')
max(max(abs(b-c)))

И если результаты (которые должны быть нулевыми) для pseudo ещё более-менее приемлемы -- в пределах семи знаков, то у pinv выходит полная катастрофа:
код: [ скачать ] [ спрятать ]
  1. B=pseudo(A) 
  2. |ABA-A|= 
  3. ans =  4.8465e-007 
  4. ans =  4.8465e-007 
  5. ans =  4.8465e-007 
  6. |BAB-B|= 
  7. ans =  2.1176e-022 
  8. ans =  4.0292e-017 
  9. ans =  3.9228e-017 
  10. |A'B'-BA|= 
  11. ans =  1.2705e-021 
  12. |B'A'-AB|= 
  13. ans =  7.1812e-013 
  14.  
  15. ================================ 
  16. C=pinv(A) 
  17. |ACA-A|= 
  18. ans =  4.8850e-015 
  19. ans =  5.0059e-015 
  20. ans =  5.0058e-015 
  21. |CAC-C|= 
  22. ans =    0.1250 
  23. ans =    0.8406 
  24. ans =    0.1366 
  25. |A'C'-CA|= 
  26. ans =  1.4901e-008 
  27. |C'A'-AC|= 
  28. ans =  1.9207e-008 
  29.  
  30. ================================ 
  31. |B|= 
  32. ans =    0.5811 
  33. |C|= 
  34. ans =  2.1429e+014 
  35. |B-C|= 
  36. ans =  2.1429e+014 


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

код: [ скачать ] [ спрятать ]
  1. B=pseudo(A) 
  2. |ABA-A|= 
  3. ans =  4.6665e-015 
  4. ans =  4.6665e-015 
  5. ans =  4.6665e-015 
  6. |BAB-B|= 
  7. ans =  3.2596e-009 
  8. ans =  3.2803e-009 
  9. ans =  3.2618e-009 
  10. |A'B'-BA|= 
  11. ans =  1.1720e-009 
  12. |B'A'-AB|= 
  13. ans =  3.7773e-015 
  14.  
  15. ================================ 
  16. C=pinv(A) 
  17. |ACA-A|= 
  18. ans =  4.6665e-015 
  19. ans =  4.6665e-015 
  20. ans =  4.6665e-015 
  21. |CAC-C|= 
  22. ans =  6.9849e-010 
  23. ans =  8.4307e-010 
  24. ans =  8.4309e-010 
  25. |A'C'-CA|= 
  26. ans =  1.9491e-015 
  27. |C'A'-AC|= 
  28. ans =  2.1266e-015 
  29.  
  30. ================================ 
  31. |B|= 
  32. ans =  2.0634e+006 
  33. |C|= 
  34. ans =  2.0634e+006 
  35. |B-C|= 
  36. ans =  7.8306e-009 

А дело действительно в плохой обусловленности матрицы:

Используется синтаксис Matlab M
a=[ 4097             -3372.3292364      5730820823.8     -14151482820       1.4429140017e+16 -5.9384620888e+16  4.3249891776e+22
   -3372.3292364      5730820823.8     -14151482820       1.4429140017e+16 -5.9384620888e+16  4.3249891776e+22 -2.4919888057e+23
    5730820823.8     -14151482820       1.4429140017e+16 -5.9384620888e+16  4.3249891776e+22 -2.4919888057e+23  1.4116057367e+29      
   -14151482820       1.4429140017e+16 -5.9384620888e+16  4.3249891776e+22 -2.4919888057e+23  1.4116057367e+29 -1.0457268159e+30
    1.4429140017e+16 -5.9384620888e+16  4.3249891776e+22 -2.4919888057e+23  1.4116057367e+29 -1.0457268159e+30  4.8465906922e+35
   -5.9384620888e+16  4.3249891776e+22 -2.4919888057e+23  1.4116057367e+29 -1.0457268159e+30  4.8465906922e+35 -4.3882409082e+36
    4.3249891776e+22 -2.4919888057e+23  1.4116057367e+29 -1.0457268159e+30  4.8465906922e+35 -4.3882409082e+36  1.7209134454e+42];

a=a*1e-42;

singeig=sort(sqrt(eig(a'*a)))
singeig=sort(sqrt(abs(eig(a'*a))))
singsvd=sort(svd(a))
err=max(abs(singeig-singsvd))



код: [ скачать ] [ спрятать ]
  1. singeig = 
  2.      1.143661295243980e-027                          
  3.                                       0 +6.494726427148507e-027i 
  4.      3.950610116881154e-021                          
  5.      2.148173051969610e-020                          
  6.                                       0 +5.329070518207753e-015i 
  7.      4.846478794280131e-007                          
  8.      1.720913445411326e+000                          
  9.  
  10. singeig = 
  11.     1.143661295243980e-027 
  12.     6.494726427148507e-027 
  13.     3.950610116881154e-021 
  14.     2.148173051969610e-020 
  15.     5.329070518207753e-015 
  16.     4.846478794280131e-007 
  17.     1.720913445411326e+000 
  18.  
  19. singsvd = 
  20.     4.988146426141484e-034 
  21.     5.845467602932113e-028 
  22.     1.553743007222084e-025 
  23.     2.135767307981534e-021 
  24.     4.666467958383161e-015 
  25.     4.846478794280069e-007 
  26.     1.720913445411327e+000 
  27.  
  28. err =    6.626025598245917e-016 

И если предпоследнее собственное число -- порядка десять в минус седьмой -- ещё следует учитывать, то третье снизу уже явно равно нулю в пределах погрешностей округления. У меня процедура по умолчанию считала ненулевым только главное сингулярное число, т.е. сочла ранг матрицы единичным (поскольку она рабодила с матрицей $A'A$ и, следовательно, неявно ориентировалось на квадраты сингулярных чисел). Вот именно с такой точностью -- до семи знаков -- она результат и выдала. Функция же pinv по наивности сочла ранг равным трём и попыталась ортогонализовать и к третьему базисному вектору тоже; а поскольку точность этого третьего вектора заведомо никакая -- то и результат вышел никакой.

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

 Профиль  
                  
 
 Re: Обращение и псевдо обращение матриц
Сообщение08.07.2012, 20:40 
Аватара пользователя


06/07/12
70
Вот сингулярные числа найденные мною с четверной точностью:
1. 0.856189752969669569039751117306343D+03
2. 0.299404523295453987913981289407457D+09
3. 0.232669668764825313446744561491917D+15
4. 0.213576730798142684509833361479408D+22
5. 0.466646795838319150077138173301492D+28
6. 0.484647879428006760004447110065820D+36
7. 0.172091344541132628614006575102588D+43

а вот сингулярные числа, найденные Вами с двойной точностью:
1. 4.988146426141484e-034
2. 5.845467602932113e-028
3. 1.553743007222084e-025
4. 2.135767307981534e-021
5. 4.666467958383161e-015
6. 4.846478794280069e-007
7. 1.720913445411327e+000

А теперь Ваше sqrt(abs(eig(a'*a))):
1. 1.143661295243980e-027
2. 6.494726427148507e-027
3. 3.950610116881154e-021
4. 2.148173051969610e-020
5. 5.329070518207753e-015
6. 4.846478794280131e-007
7. 1.720913445411326e+000

Т.е. если 5-ые сингулярные числа совпадают с точностью до14 знака (интересный вопрос: почему? Ответ кроется в структуре матрицы: сначала обрабатываются небольшие (относительно!) элементы, и так по нарастающей),
то 5-ое singeig уже никуда не годится. Почему? Ответ очевиден: перемножаете 2-е плохие м-цы и получаете очень плохую матрицу. И потом я не понял Вашу фразу: "попыталась ортогонализовать".

 Профиль  
                  
 
 Re: Обращение и псевдо обращение матриц
Сообщение08.07.2012, 22:15 
Заслуженный участник


11/05/08
31922
Да тут бессмысленно применять четверную точность, поскольку исходные данные в грубейшем элементе заданы лишь четырьмя значащими цифрами, да и остальные -- максимум одиннадцатью, что многократно грубее даже точности double. Собственно, даже учёт второго сингулярного числа, приводящий к несколько экзотическому результату -- выглядит в этой ситуации весьма сомнительным. Слишком уж велика погрешность входных данных, чтобы можно было надеяться на хоть мало-мальскую устойчивость результата при учёте младших сингулярных чисел. Т.е. результат-то мы получим, конечно; но показывать он будет среднюю температуру по больнице.

-- Вс июл 08, 2012 23:48:07 --

Добивка. Почему я считаю учёт уже второго сингулярного числа в этой ситуации сомнительным?

Во-первых, как уже сказал -- потому, что он приводит к слишком уж сильному (на несколько порядков) изменению результата. Но не только поэтому.

Я ведь сознательно привёл сингулярные числа, вычисленные двумя способами: и стандартной svd, и тупо через собственные числа $A'A$. Главные сингулярные числа при этом, естественно, совпали (кто бы и сомневался). А вот то, что практически совпали и вторые (будучи при этом маленькими) -- это уже наводит на определённые размышления.

Это наводит на мысль о том, что процедура svd работает фактически тоже с матрицей $A'A$. Между тем для последней разброс между первым и вторым собственными числами составляет уже 13 порядков, что многократно меньше относительной погрешности входных данных. Откуда и вывод: второе сингулярное число для данной матрицы -- это наверняка филькина грамота. Не говоря уж о следующих.

 Профиль  
                  
 
 Re: Обращение и псевдо обращение матриц
Сообщение08.07.2012, 23:03 
Аватара пользователя


06/07/12
70
Откиньте сколько Вам угодно знаков - вывод не изменится: SVD (правда, в данном конкретном случае достаточно и диагонализации) не переплюнуть (вообще говоря, я не нашел в тексте место, в котором говорится о точности задания элементов матрицы). И по поводу 4-ой точности в данном конкретном случае Вы не правы. Четверная точность здесь полезна тем, что сравнивая решения с двойной и четверной точностью, мы можем судить о числе правильных значащих цифр, найденных с 2-ой точностью: распределение элементов в матрице таково, что найденные с.в. и с.з. находятся со значительно более высокой точностью, чем ожидается на первый взгляд и четверная точность сильно облегчает анализ этого феномена. Кстати, многие трехдиагональные симметричные матрицы в мат. физике имеют аналогичную структуру с монотонным возрастанием элементов, при этом относительный разброс может достигать громадных значений, а с.в. и с.з. при этом будут найдены очень точно.

-- 08.07.2012, 23:12 --

ewert в сообщении #593630 писал(а):
А вот то, что практически совпали и вторые (будучи при этом маленькими) -- это уже наводит на определённые размышления.

Это наводит на мысль о том, что процедура svd работает фактически тоже с матрицей $A'A$. Между тем для последней разброс между первым и вторым собственными числами составляет уже 13 порядков, что многократно меньше относительной погрешности входных данных. Откуда и вывод: второе сингулярное число для данной матрицы -- это наверняка филькина грамота. Не говоря уж о следующих.

Вы не правы и выше объяснено, почему. А про $A'A$ - без комментариев.

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

Модераторы: Toucan, maxal, PAV, Karan, Супермодераторы



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

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


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

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