АЛГОРИТМ ХАУСХОЛДЕРА

Преобразование Хаусхолдера. Матрица преобразования Хаусхолдера описывается резольвентой Q=(E-2uu'/u'u).

Это преобразование инвертирует вектор u в том смысле, что Qu=-u, а любой ортогональный к нему вектор s не затрагивает Qs=s. В трехмерном пространстве это зеркальное отражение вектора u. По отношению к произвольному вектору v задача ставится так, чтобы найти u, приводящее его к масштабированному орту e, то есть Qv=ne или Qv=-ne, все равно, где n=√(v'v). Знак при норме связывают со значением отвечающего единичному элементу e элементу v, для отрицательных элементов этот знак положительный.

Алгоритм Хаусхолдера. Более общая задача состоит в том, чтобы анулировать только часть нижних компонент вектора v, начиная с индекса p, меняя при этом главный элемент (ведущий) под номером j. В таком случае норма n вычисляется с усеченным вектором, включающим только изменяемые компоненты v. Знак ведущего элемента меняется на противоположный. Задача оперирует относительно небольшим количеством переменных, поэтому алгоритм Хаусхолдера состоит в применении установленных им формул расчета. Применим эти формулы и покажем, что они приводят к нужному результату.


v=rand(10), j=2, p=6, 
u=H(j p v), Q=getQ(u), y={Q*v}, mesh(Q), ball(y),

function: H(j p v),
var i s, j=j-1, p=p-1, 
u={v}, for i=0:p-1, u[i]=0, end, 
s=v[j]*v[j], for i=p:size(v), s=s+v[i]*v[i], end, 
s=sqrt(s), if v[j]>0, s=-s, end, u[j]=v[j]-s, 
return u,
end,

function: getQ(u),
I=eye(rows(u)), u2={u'*u}, u2={u/u2}, u2={-2*u2}, 
Q=res(I+u2*u), 
return Q,
end,

Модифицированный алгоритм Хаусхолдера. Матрица Q преобразования Хаусхолдера, экономя память, в явном виде не вычисляют, а вычисляют сразу итог ее умножения на матрицу С, сохраняя результат в той же матрице (если начинать с единичной матрицы, то на ее месте получим Q). Вектор u содержит ведущий элемент, нулевые элементы и элементы исходного вектора v, начиная с номера k. Как видно, достаточно рассчитывать только ведущий элемент в u. Вектор y=Qv наоборот, отличен от нуля только в верхней его части, которую сохраняют в v.


j=2, p=6, v=rand(10), C=eye(10), v2={v}, C2={C}, 

call u=H2(j p v C), mesh(C), u=?,
% СРАВНЕНИЕ С ПОЭТАПНЫМ ВЫЧИСЛЕНИЕМ,
call u=H(j p v2), Q=getQ(u), 
C3={Q*C2}, n={C3-C}, n=norm(n), n=?,

function: H2(j p v C),
var i k s f, j=j-1, p=p-1,  
s=v[j]*v[j], for i=p:size(v), s=s+v[i]*v[i], end, 
s=sqrt(s), if v[j]>0, s=-s, end, u=v[j]-s, v[j]=s,
f=s*u, if f<>0, 
for i=0:colm(C), s=C[i][j]*u, 
for k=p:rowm(C), s=s+C[i][k]*v[k], end, s=s/f,
C[i][j]=C[i][j]+s*u, 
for k=p:rowm(C), C[i][k]=C[i][k]+s*v[k], end,
end, end, 
return u,
end,

function: H(j p v),
var i s, j=j-1, p=p-1, 
u={v}, for i=0:p-1, u[i]=0, end, 
s=v[j]*v[j], for i=p:size(v), s=s+v[i]*v[i], end, 
s=sqrt(s), if v[j]>0, s=-s, end, u[j]=v[j]-s, 
return u,
end,

function: getQ(u),
I=eye(rows(u)), u2={u'*u}, u2={u/u2}, u2={-2*u2}, 
Q=res(I+u2*u), 
return Q,
end,

На основе этого алгоритма пишется QR разложнение матрицы P=Q'R, где Q представляет собой ортогональную матрицу, а R - это верхняя треугольная матрица.

Содержание алгоритма довольно легко пояснить, опираясь на свойства преобразования Хаусхолдера. Сначала к единичному орту преобразуется первый столбец матрицы P. Остальные столбцы преобразуются так, как преобразуются. Повлиять на них направленно мы уже не можем. Потом применяется блочно составная матрица преобразования Хаусхолдера, начало которой начинает от преобразования к преобразованию замещаться растущей единичной матрицей. Благодаря этому верхний элемент полученного орта остается без изменений, анулированная часть преобразованием не меняется, то есть, первый столбец таким преобразованием сохраняется. Преобразование нацеливается на то, чтобы анулировать друг за другом низ все укорачивающихся столбцов. Итогом будет верхняя треугольная матрица R.



Rambler's Top100