|
|
Ну пожалуйста. Сначала программируем нахождение псевдорешения (какого-либо) системы  :
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 к ядру:
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
Вот тестовая программка:
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. Правда, точность на пару знаков ниже, но я ведь и не пытался выжать все соки. Тем более не пытался что-то оптимизировать.
|
|
04.07.2012, 12:03 
правых частей нормальное псевдорешение, т.е. решение 
системы 
, имеющее минимальную норму. Собственно решение можно найти, скажем, обычным методом Гаусса. Причём если решений много, то достаточно найти одно; тогда нормальное решение получается просто ортогонализацией этого одного решения ядру матрицы 
.
Что такое 
называется решение, имеющее наименьшую норму (если, конечно, решений бесконечно много). Или, что эквивалентно -- решение, ортогональное ядру матрицы.
называется матрица 
, которая каждый вектор 
, являющийся нормальным псевдорешением системы 
и, следовательно, неявно ориентировалось на