Prvni kroky bidiagonalizace nasledujici matice:

>> A = [3, -5, 15, -4; 7, 7, 9 -1; 0, 8, -1, 6; -10, 7, 8, -10]
A =

    3   -5   15   -4
    7    7    9   -1
    0    8   -1    6
  -10    7    8  -10



Vezmu prvni sloupec, potrebuji snulovat vsechny slozky krome prvni.

>> x = A(:,1)
x =

    3
    7
    0
  -10

Vytvorim, stejne velky vektor, ktery ma ty prvky nulove:

>>  y = [norm(x);0;0;0]
y =

   12.56981
    0.00000
    0.00000
    0.00000


Vypocitam Householderovu transformaci:

>> w = (x-y)/norm(x-y)
w =

  -0.61698
   0.45130
   0.00000
  -0.64472

>> H = eye(4) - 2*w*w'
H =

   0.23867   0.55689  -0.00000  -0.79556
   0.55689   0.59265  -0.00000   0.58192
  -0.00000  -0.00000   1.00000  -0.00000
  -0.79556   0.58192  -0.00000   0.16868


Vypoctena Householderova transformace nuluje mimo(bi)diagonalni prvky v prvnim sloupci:

>> H*A
ans =

   12.56981   -2.86401    2.22756    6.44401
    0.00000    5.43759   18.34262   -8.63946
    0.00000    8.00000   -1.00000    6.00000
    0.00000    9.23201   -5.34660    0.91351

Prvni sloupec je tedy ok, odrezu ho, a zabyvam se transponovanou matici:

>> A1 = (H*A)'(2:4,:)
A1 =

   -2.86401    5.43759    8.00000    9.23201
    2.22756   18.34262   -1.00000   -5.34660
    6.44401   -8.63946    6.00000    0.91351


Stejny postup: vezmu prvni sloupec a stejne velky vektor se snulovanymi slozkami (mimo prvni).

>> x = A1(:,1)
x =

  -2.8640
   2.2276
   6.4440


>>  y = [norm(x);0;0]
y =

   7.39526
   0.00000
   0.00000

Vypoctu Householderovu transformaci:

>> w = (x-y)/norm(x-y)
w =

  -0.83285
   0.18083
   0.52313

>> W = eye(3) - 2*w*w'
W =

  -0.38728   0.30121   0.87137
   0.30121   0.93460  -0.18920
   0.87137  -0.18920   0.45268


Tato transformace nuluje mimo(bi)diagonalni prvky v prvnim sloupci
(nevychazi to presne kvuli zaokrouhlovaci chybe).

>> W * A1
ans =

  7.3953e+000  -4.1090e+000  1.8288e+000  -4.3898e+000
  -2.2204e-016  2.0415e+001  3.3993e-001  -2.3889e+000
  -4.4409e-016  -2.6431e+000  9.8762e+000  9.4696e+000


Upravime vypoctenou transforamaci tak, aby se aplikovala nasobenim zprava na puvodni matici:

>> WT = prepad(prepad(W,4)',4);
>> WT(1,1)=1
WT =

   1.00000   0.00000   0.00000   0.00000
   0.00000  -0.38728   0.30121   0.87137
   0.00000   0.30121   0.93460  -0.18920
   0.00000   0.87137  -0.18920   0.45268

Mame poreseny prvni sloupec a prvni radek:

>> H * A * WT
ans =

   12.56981    7.39526   -0.00000   -0.00000
    0.00000   -4.10895   20.41543   -2.64312
    0.00000    1.82880    0.33993    9.87623
    0.00000   -4.38981   -2.38894    9.46960

Pokracujeme stejnym zpusobem pro submatici:

>> A1 = (H * A * WT)(2:4,2:4)
A1 =

   -4.10895   20.41543   -2.64312
    1.82880    0.33993    9.87623
   -4.38981   -2.38894    9.46960