QR-разложение. Опираясь на предыдущие алгоритмы, предложим вариант реализации ортогонального разложения матрицы, у которой число столбцов равно или меньше числа строк - со строчным преобладанием: R=QP, P=Q'R, причем R формируется по месту P.


P=[1 2; 3 4; 5 6], P2={P}, 

Q=qr2(P), mesh(P), 
% ПРОВЕРКА ИТОГОВ РАЗЛОЖЕНИЯ,
P2={Q*P2}, P2={P2-P}, dif=norm(P2), dif=?,     

function: qr2(R),
var i j k n m u v f,
m=rowm(R), n=colm(R), Q=eye(m+1),      
for j=0:n, 
  s=0.0, for i=j:m, v=R[j][i], s=s+v*v, end, s=sqrt(s),
  if R[j][j]>0, s=-s, end, u=R[j][j]-s, R[j][j]=s,
  f=R[j][j]*u,
  if f<>0, for i=j+1:n, s=R[i][j]*u,
           for k=j+1:m, s=s+R[i][k]*R[j][k], end, s=s/f,
           R[i][j]=R[i][j]+s*u,
           for k=j+1:m, R[i][k]=R[i][k]+s*R[j][k], end,
           end,
           for i=0:m, s=Q[i][j]*u,
           for k=j+1:m, s=s+Q[i][k]*R[j][k], end, s=s/f,           
           Q[i][j]=Q[i][j]+s*u,
           for k=j+1:m, Q[i][k]=Q[i][k]+s*R[j][k], end,
           end,
  end,  
end,
  for i=0:m, for j=0:i-1, if j<=n, R[j][i]=0.0, end, end, end,    
return Q,
end,

МЕТОД НАИМЕНЬШИХ КВАДРАТОВ

Задача решения системы Ax=b методом наименьших квадратов приводит к системе нормальных уравнений Px=r, где P=A'A, r=A'b. Невырожденные уравнения можно решать, используя разложение Холецкого P=LL'. В более общем случае возможно выполнить трапецевидное разложение P=TT', переставляя анулированные столбцы в конец.

Обратим внимание, что транспонированная к T матрица в явном виде показывает уравнения, которые можно отбраковать.

Усеченная система имеет вид T'x=L-1r, где L в данном случае - треугольник на вершине трапеции, T - прямоугольная матрица, полученная отбрасыванием нулевых столбцов. Модифицированный алгоритм Холецкого возвращает, как и в случае полного ранга, правую часть в готовом виде. Поэтому наши заботы сводятся теперь к триангуляризации прямоугольной матрицы T'.

Запишем усеченную систему в виде T'Q'Qx=L-1r, где Q - матрица преобразования Хаусхолдера, приводящая T к верхнему треугольному виду R=QT. Правая трапецевидная матрица R отличается от левой трапеции T тем, что нижние строки ее анулированы. Уравнение сводится к R'Qx=L-1r.

Поступим также, как и раньше, отбрасывая нулевую кайму в R, тогда Q сократится до прямоугольного блока с тем же самым количеством строчек, для согласования размерностей. Матрица R становится попросту треугольной, наконец, мы избавились от трапеций. Решение x=Q'R'-1L-1r.


A=rands(4), b=[1 1 1 1], a=rand(4), A[0]=a, A[1]=a,
% ГОТОВИМ ВЫРОЖДЕННУЮ МАТРИЦУ,  
P={A'*A}, r={A'*b}, I=line(r), P2={P}, r2={r}, 
% ТРАПЕЦЕВИДНОЕ РАЗЛОЖЕНИЕ, 
call rank=ch5(P r I 0.001), call P=perm(P r I rank),
% ГРАФИК ПАДЕНИЯ НЕВЯЗКИ МНК,    
n=norm(b), n=n*n, g=one(rank),    
for i=0:rank-1, v=r[i], n=abs(n-v*v), g[i]=n, end, ball(g), rank=?,
% ПРЕОБРАЗОВАНИЕ ХАУСХОЛДЕРА и УСЕЧЕНИЯ,
Q=qr2(P), 
L=cutrR(P rank), mesh(L),
Q=cutrQ(Q rank), 
% РЕШЕНИЕ ТРЕУГОЛЬНОЙ СИСТЕМЫ,
x=triangle(L r), x={Q*x},
% ОБРАТНАЯ ПЕРЕСТАНОВКА, 
x=restorex(x I), % x=?,  

% ПРОВЕРКА РЕЗУЛЬТАТА С БИБЛИОТЕЧНОЙ ФУНКЦИЕЙ,
xx={A\b}, dif={x-xx}, dif=norm(dif), dif=?, 

function: ch5(L r I eps),
var ii jj kk i j k s max,
s=rowm(L), rank=s+1,  
max=reset(L r I 0 eps), 
if max<eps*eps, rank=0, else, kk=I[0],
L[kk][kk]=max, r[kk]=r[kk]/max, 
for i=1:s, ii=I[i], L[kk][ii]=L[kk][ii]/max, L[ii][kk]=0, 
L[ii][ii]=L[ii][ii]-L[kk][ii]*L[kk][ii],  
r[ii]=r[ii]-L[kk][ii]*r[kk],
end,
for k=1:s, if rank=s+1,  
max=reset(L r I k eps), 
if max<eps*eps, rank=k, else, kk=I[k], L[kk][kk]=max,
for i=k+1:s, ii=I[i], for j=0:k-1, jj=I[j], 
L[kk][ii]=L[kk][ii]-L[jj][ii]*L[jj][kk], end, 
L[kk][ii]=L[kk][ii]/max, end, 
r[kk]=r[kk]/max,
for i=k+1:s, ii=I[i],
L[ii][ii]=L[ii][ii]-L[kk][ii]*L[kk][ii], L[ii][kk]=0,
r[ii]=r[ii]-L[kk][ii]*r[kk], 
end, end, end, end, end,
return rank,
end,

function: reset(L r I k eps), 
var eps imax rr Lii ii i s,
s=size(I), imax=k, max=0, 
% ВЫБИРАЕТСЯ МАКСИМАЛЬНОЕ ПАДЕНИЕ НЕВЯЗКИ,
for i=k:s, ii=I[i], rr=abs(r[ii]), Lii=sqrt(abs(L[ii][ii])), 
if Lii>eps*rr, if rr>max*Lii, max=rr/Lii, imax=i, end, end, end,
if max=0, imax=k,  
% ВЫБИРАЕТСЯ МАКСИМАЛЬНЫЙ ДИАГОНАЛЬНЫЙ ЭЛЕМЕНТ,
for i=k:s, ii=I[i], Lii=sqrt(abs(L[ii][ii])), 
if Lii>max, max=Lii, imax=i, end, end,  
end,
ii=I[imax], I[imax]=I[k], I[k]=ii, max=sqrt(abs(L[ii][ii])), 
return max,
end,

function: perm(L r I rank),
var s i ii buf, s=size(I), 
buf={r}, for i=0:s, ii=I[i], r[i]=buf[ii], end, 
% ОТРЕЗАЕМ ЛИШНИЕ ЭЛЕМЕНТЫ ХВОСТА,
for i=0:s-rank, r.pop(), end,
buf={L}, L=one(1), for i=0:rank-1, ii=I[i], L[i]={buf[ii]}, end,
buf={L}, tr(buf),  for i=0:s, ii=I[i], L[i]={buf[ii]}, end, 
tr(L), return L,
end,

function: qr2(R),
var i j k n m u v f,
m=rows(R), n=cols(R), Q=eye(m), n=n-1, m=m-1,    
for j=0:n, 
  s=0.0, for i=j:m, v=R[j][i], s=s+v*v, end, s=sqrt(s),
  if R[j][j]>0, s=-s, end, u=R[j][j]-s, R[j][j]=s,
  f=R[j][j]*u,
  if f<>0, for i=j+1:n, s=R[i][j]*u,
           for k=j+1:m, s=s+R[i][k]*R[j][k], end, s=s/f,
           R[i][j]=R[i][j]+s*u,
           for k=j+1:m, R[i][k]=R[i][k]+s*R[j][k], end,
           end,
           for i=0:m, s=Q[i][j]*u,
           for k=j+1:m, s=s+Q[i][k]*R[j][k], end, s=s/f,           
           Q[i][j]=Q[i][j]+s*u,
           for k=j+1:m, Q[i][k]=Q[i][k]+s*R[j][k], end,
           end,
  end,  
end,
  for i=0:m, for j=0:i-1, if j<=n, R[j][i]=0.0, end, end, end,    
return Q,
end,

function: cutrR(R rank),
var i j buf,   
buf={R}, R=ones(rank),  
for i=0:rank-1, for j=0:rank-1, R[i][j]=buf[j][i], end, end,
return R,
end,

function: cutrQ(Q rank),
var i j buf,   
buf={Q}, Q=ones(rank),  
for i=0:rank-1, for j=0:size(buf), Q[i][j]=buf[j][i], end, end,  
return Q,
end,

function: triangle(L r),
var s i j, x=one(r), x[0]=r[0]/L[0][0],
for i=1:size(x), s=r[i], 
for j=0:i-1, s=s-L[j][i]*x[j], end, x[i]=s/L[i][i], end,  
return x,
end,

function: restorex(x I),
var s i ii buf, s=size(I), buf={x}, 
for i=0:s, ii=I[i], x[ii]=buf[i], end,  
return x,
end,

В завершение учтен индексный массив. Остается собрать составные части, чтобы посмотреть примеры.



Rambler's Top100