АЛГОРИТМЫ
Новости
Рассылка новостей
Форум
AlgoPascal
Редактор блок-схем
Статьи
О сайте
Контакты
Содержание - Операции с матрицами и векторами - Сингулярное разложение (SVD-разложение) произвольной матрицыСингулярное разложение (SVD-разложение) произвольной матрицы
Сингулярное разложение матрицы (ещё его называют SVD-разложением) - это представление произвольной матрицы A (размер M*N) в виде произведения прямоугольной матрицы U (размер M*N, с ортонормированными или нулевыми столбцами), диагональной матрицы W (размер N*N) с неотрицательными элементами на главной диагонали и квадратной матрицы V (размер N*N) с ортонормированными столбцами:
A = U*diag(w1 , ..., wn )*V TПредставление матрицы в такой форме позволяет решать целый ряд задач. Кратко перечислю некоторые полезные свойства SVD-разложения:
- Число обусловленности матрицы равно отношению максимального из wi и минимального.
- Столбцы матрицы U, чьим номерам соответсвуют ненулевые элементы wi , образуют ортонормированный базис линейного подпространства, образованного столбцами матрицы A.
- Столбцы матрицы V, чьим номерам соответсвуют нулевые элементы wi , образуют ортонормированный базис линейного подпространства ненулевых решений системы уравнений Ax = 0.
- Матрица, обратная к матрице A, имеет вид: A -1 = V*diag(1/wi )*U T
- Если матрица A вырожденная, то наилучшее приближение к решению системы линейных уравнений Ax = b ищется по предыдущей формуле, с умножением справа на b, причем те элементы матрицы diag(1/wi ), в знаменателе которых стоит ноль, мы полагаем равными нолю (это не ошибка!).
Алгоритм, осуществляющий сингулярное разложение A, имеет трудоемкость порядка N 3. Он является итерационным, но даже на матрицах порядка 500*500 обычно требуется не более 8 итераций. Ходят слухи, что в некоторых редчайших случаях этот алгоритм не сходится за отведенные с запасом несколько десятков итераций. В таком случае функция возвращает False, но такое происходит очень редко и я пока ни разу не получал такой эффект. В общем, это не то, из-за чего следует беспокоиться. Алгоритм очень устойчив, и вы вряд ли когда-нибудь получите результат, свидетельствующий о несходимости итераций.
Если нашли ошибку в алгоритме - сообщите!
Реализации алгоритма на различных языках:
Реализация алгоритма на C++Реализация алгоритма на DelphiРеализация алгоритма на Visual Basic 6Реализация алгоритма на AlgoPascal:
unit SVDDecompositionUnit;
interface SVDDecomposition;
implementation
const
(***********************************************
Максимальное число внутренних итераций алгоритма
Обычно сходимость уже на 7-9 итерациях. Взято с
запасом, благо скорость от этого не падает.
***********************************************)
MaxSVDIterations : Integer = 60;
(***********************************************************************
Алгоритм сингулярного разложения матрицы размером MxN.
Алгоритм принимает матрицу A и приводит её к виду: A = U*W*Transpone(V).
Матрица U имеет размер MxN, матрица W - диагональная, матрица V орто-
нормированная.
Входные данные:
*Матрица А, строки с 1 по M, столбцы с 1 по N.
*M, N - размер матрицы.
Результат работы помещается:
*Матрица U замещает матрицу A, строки с 1 по M, столбцы с 1 по N.
*Диагональ матрицы W хранится в переменной W (элементы с 1 по N).
*Матрица V (не транспонированная) хранится в переменной V. Строки
с 1 по N, столбцы с 1 по N.
Если алгоритм закончил работу корректно, то функция возвращает True. Если
же итерации алгоритма не привели к сходимости (редчайший случай), то воз
вращается False;
***********************************************************************)
function SVDDecomposition(
var a : array of array of Real;
m : Integer;
n : Integer;
out w : array of Real;
out v : array of array of Real):Boolean;
var
nm : Integer;
MinMN: Integer;
l : Integer;
k : Integer;
j : Integer;
jj : Integer;
its : Integer;
i : Integer;
z : Real;
y : Real;
x : Real;
scale: Real;
s : Real;
h : Real;
g : Real;
f : Real;
c : Real;
anorm: Real;
rv1: array of Real;
Flag: Boolean;
begin
SetBounds(rv1, [1, N]);
SetBounds(W, [1,N]);
SetBounds(V, [1,N], [1,N]);
Result:=True;
if M<N then
MinMN:=M
else
MinMN:=N;
g := 0.0;
scale:=0.0;
anorm:=0.0;
for i:=1 to n do
begin
l:=i+1;
rv1[i]:=scale*g;
g:=0;
s:=0;
scale:=0;
if i<=m then
begin
for k:=i to m do
scale:=scale+AbsReal(a[k,i]);
if scale<>0.0 then
begin
for k:=i to m do
begin
a[k,i] := a[k,i]/scale;
s := s+a[k,i]*a[k,i]
end;
f:=a[i,i];
g:=-ExtSign(sqrt(s),f);
h:=f*g-s;
a[i,i]:=f-g;
if i<>n then
begin
for j:=l to n do
begin
s:=0.0;
for k:=i to m do
s:=s+a[k,i]*a[k,j];
f:=s/h;
for k:=i to m do
a[k,j] := a[k,j]+f*a[k,i];
end;
end;
for k:=i to m do
a[k,i] := scale*a[k,i];
end;
end;
w[i]:=scale*g;
g:=0.0;
s:=0.0;
scale:=0.0;
if (i<=m) and (i<>n) then
begin
for k:=l to n do
scale:=scale+AbsReal(a[i,k]);
if scale<>0.0 then
begin
for k:=l to n do
begin
a[i,k]:=a[i,k]/scale;
s:=s+a[i,k]*a[i,k];
end;
f:=a[i,l];
g:=-ExtSign(sqrt(s),f);
h:=f*g-s;
a[i,l]:=f-g;
for k:=l to n do
rv1[k] := a[i,k]/h;
if i<>m then
begin
for j:=l to m do
begin
s:=0.0;
for k:=l to n do
s:=s+a[j,k]*a[i,k];
for k:=l to n do
a[j,k]:=a[j,k]+s*rv1[k];
end;
end;
for k:=l to n do
a[i,k]:=scale*a[i,k];
end;
end;
anorm:=MyMax(anorm,(AbsReal(w[i])+AbsReal(rv1[i])))
end;
for i:=n downto 1 do
begin
if i<n then
begin
if g<>0.0 then
begin
for j:=l to n do
v[j,i] := (a[i,j]/a[i,l])/g;
for j:=l to n do
begin
s:=0.0;
for k:=l to n do
s:=s+a[i,k]*v[k,j];
for k:=l to n do
v[k,j]:=v[k,j]+s*v[k,i];
end;
end;
for j:=l to n do
begin
v[i,j]:=0.0;
v[j,i]:=0.0;
end;
end;
v[i,i]:=1.0;
g:=rv1[i];
l:=i;
end;
for i:=MinMN downto 1 do
begin
l:=i+1;
g:=w[i];
if i<n then
begin
for j:=l to n do
a[i,j]:=0.0;
end;
if g<>0.0 then
begin
g:=1.0/g;
if i<>n then
begin
for j:=l to n do
begin
s := 0.0;
for k:=l to m do
s:=s+a[k,i]*a[k,j];
f:=(s/a[i,i])*g;
for k:=i to m do
a[k,j]:=a[k,j]+f*a[k,i];
end;
end;
for j:=i to m do
a[j,i]:=a[j,i]*g;
end
else
begin
for j:=i to m do
a[j,i]:=0.0;
end;
a[i,i]:=a[i,i]+1.0;
end;
for k:=n downto 1 do
begin
for its:=1 to MaxSVDIterations do
begin
Flag:=True;
for l:=k downto 1 do
begin
nm:=l-1;
if (AbsReal(rv1[l])+anorm)=anorm then
begin
Flag:=False;
Break;
end;
if (AbsReal(w[nm])+anorm)=anorm then
begin
Break;
end;
end;
if Flag then
begin
c:=0.0;
s:=1.0;
for i:=l to k do
begin
f:=s*rv1[i];
if (AbsReal(f)+anorm)<>anorm then
begin
g:=w[i];
h:=Pythag(f,g);
w[i]:=h;
h:=1.0/h;
c:=g*h;
s:=-(f*h);
for j:=1 to m do
begin
y:=a[j,nm];
z:=a[j,i];
a[j,nm]:=(y*c)+(z*s);
a[j,i]:=-(y*s)+(z*c);
end;
end;
end;
end;
z:=w[k];
if l=k then
begin
if z<0.0 then
begin
w[k]:=-z;
for j:=1 to n do
v[j,k]:=-v[j,k];
end;
Break;
end;
if its=MaxSVDIterations then
begin
Result:=False;
Exit;
end;
x:=w[l];
nm:=k-1;
y:=w[nm];
g:=rv1[nm];
h:=rv1[k];
f:=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
g:=Pythag(f,1);
f:=((x-z)*(x+z)+h*((y/(f+ExtSign(g,f)))-h))/x;
c:=1.0;
s:=1.0;
for j:=l to nm do
begin
i:=j+1;
g:=rv1[i];
y:=w[i];
h:=s*g;
g:=c*g;
z:=Pythag(f,h);
rv1[j]:=z;
c:=f/z;
s:=h/z;
f:=(x*c)+(g*s);
g:=-(x*s)+(g*c);
h:=y*s;
y:=y*c;
for jj:=1 to n do
begin
x:=v[jj,j];
z:=v[jj,i];
v[jj,j]:=(x*c)+(z*s);
v[jj,i]:=-(x*s)+(z*c)
end;
z:=Pythag(f,h);
w[j]:=z;
if z<>0.0 then
begin
z:=1.0/z;
c:=f*z;
s:=h*z;
end;
f:=(c*g)+(s*y);
x:=-(s*g)+(c*y);
for jj:=1 to m do
begin
y:=a[jj,j];
z:=a[jj,i];
a[jj,j]:=(y*c)+(z*s);
a[jj,i]:=-(y*s)+(z*c)
end;
end;
rv1[l]:=0.0;
rv1[k]:=f;
w[k]:=x;
end;
end;
end;
(************************************************************
Служебные подпрограммы
************************************************************)
function ExtSign(a:Real; b: Real):Real;
begin
if b>=0 then
Result:=AbsReal(a)
else
Result:=-AbsReal(a);
end;
function MyMax(a:Real; b:Real):Real;
begin
if a>b then
Result:=a
else
Result:=b;
end;
function Pythag(A:Real; B:Real):Real;
begin
if AbsReal(A)<AbsReal(B) then
Result:=AbsReal(B)*Sqrt(1+Sqr(A/B))
else
Result:=AbsReal(A)*Sqrt(1+Sqr(B/A));
end;
end.
Бочканов Сергей, Быстрицкий Владимир
Copyright © 1999-2004
При поддержке проекта MANUAL.RU