Метод Хаусхолдера, который будет описан ниже, является таким же устойчивым, как и редукция Гивенса, но требует в 2 раза меньше операций, таким образом, метод Гивенса обычно не используется.  Однако недавние исследования [6] показали, что редукция Гивенса может быть переформулирована так, чтобы повысить ее эффективность примерно в 2 раза, а также избежать вычисления квадратных корней. В результате этого "быстрый метод Гивенса" становится конкурентоспособным с методом Хаусхолдера. С другой стороны, быстрый метод должен отслеживать возможности переполнения, и переменные должны быть периодически перемасштабируемы при его использовании. Это обстоятельство говорит о том, что метод Хаусхолдера все же предпочтительнее.

Перепишем P как

Для приведения симметричной матрицы к трехдиагональной форме, вектор x, который будет образовывать первую матрицу Хаусхолдера, составим из нижних n - 1 элементов первого столбца. Тогда будут обнулены нижние n - 2 элемента:
 

Здесь матрицы записаны в блочной форме, причем ^(n-1)P₁ обозначает матрицу Хаусхолдера размерами ((n-1) x (n-1)). Величина k это просто плюс или минус модуль вектора [a₂₁, ..., a_n1]^T.

Полная ортогональная трансформация теперь будет выглядеть как:
 

Здесь был использован факт, что P^T = P.

Для второй трансформации матрица Хаусхолдера составляется на основе нижних (n - 2) элементов второго столбца, и из нее конструируется матрица преобразования:
 

Единичный блок в левом верхнем углу обеспечивает сохранение трехдиагональности, полученной на предыдущей стадии, в то время как матрица Хаусхолдера ^(n-2)P₂ размера (n-2) добавляет еще один столбец и строку к трехдиагональному результату. Ясно, что цепочка (n-2) таких преобразований приведет матрицу к трехдиагональному виду.

Вместо того, чтобы явно проводить матричные умножения в выражениях PAP, мы вычислим вектор

Программа редукции Хаусхолдера, основанная на алгоритме из [2] и помещенная ниже, реально начинает действие с n-ого столбца матрицы , а не с первого, как было написано выше. В деталях алгоритм выглядит следующим образом. На стадии m (m = 1, 2, ..., n-2) вектор u имеет вид:

Переменные вычисляются в следующем порядке: s, u, H, p, K, q, A'. На каждой m-ой стадии матрица становится трехдиагональной в последних m - 1 строках и столбцах.

Если уже найдены собственные вектора трехдиагональной матрицы (например, по программе из следующей секции), то собственные вектора могут быть найдены умножением накопленного произведения трансформирующих матриц

На входе программы задается действительная симметричная матрица в массиве a[1...n][1...n]. На выходе этот массив содержит элементы матрицы Q. Вектор d[1...n] содержит множество диагональных элементов A', вектор e[1...n] содержит внедиагональные элементы A', начиная с e[2]. Заметим, что поскольку исходное содержимое массива a уничтожается, то при производстве серии последовательных вычислений с он должен быть скопирован перед передачей в программу.

Промежуточные результаты не требуют дополнительной памяти. На стадии m вектора p и q не равны нулю только для элементов с индексами 1 ... i (вспомним, что i = n - m + 1), а вектор u -- только для элементов с индексами 1 ... i - 1. Элементы вектора e определяются в порядке n, n - 1, ..., так что вектор p можно держать на месте неопределенных элементов e. Вектор q может записываться поверх p, так как после определения q вектор p уже не требуется. Вектор u размещается в i-ой строке массива a, а u/H -- в i-ом столбце. После проведения редукции, матрицы Q_j вычисляются с использованием записанных в массив a векторов u и u/H. Поскольку матрица Q_j является единичной для последних n - j + 1 строк и столбцов, то требуется вычислить ее элементы лишь до строки и столбца n - j. При этом можно перезаписать u и u/H в соответствующих строке и столбце, поскольку они более не нужны.

Расположенная ниже программа tred2 содержит еще одну тонкость. Если величина s² равна нулю или "мала" на какой-либо из стадий, то соответствующее преобразование может быть пропущено. Простой критерий типа

Заметим, что если мы имеем дело с матрицами, элементы которых различаются на много порядков по величине, важно, чтобы строки или столбцы матрицы были переставлены, насколько это возможно, так, чтобы наименьшие элементы оказались в верхнем левом углу. Это связано с тем, что процесс начинается от правого нижнего угла, а совместные расчеты с малыми и большими величинами могут приводить к значительным ошибкам округления.

Программа tred2 предназначена для совместного использования с программой tqli, помещенной в следующей секции. tqli находит собственные вектора и значения симметричной трехдиагональной матрицы. Комбинация tred2 и tqli на сегодня является наиболее эффективным алгоритмом поиска всех собственных чисел (и, при необходимости, собственных векторов) действительной симметричной матрицы.

В листинге программы, расположенной ниже, имеются комментарии о том, какие инструкции требуются только для вычисления собственных векторов. Если закомментировать эти инструкции, то программа будет вычислять одни собственные значения, что ускорит выполнение ее приблизительно в 2 раза. Для больших значений n число операций для выполнения редукции Хаусхолдера имеет порядок 2n³/3 для вычисления одних только собственных значений, и 4n³/3 -- для собственных значений и векторов.

Программа tred2

void tred2(float **a, int n, float *d, float *e) {
  int l,k,j,i;
  float scale,hh,h,g,f;
  /* Проход по стадиям процесса редукции */
  for(i=n;i>=2;i--) {
    l=i-1; h=scale=0.;
    /* сложный процесс везде, кроме последней стадии */
    if(l>1) {
      /* вычислить шкалу */
      for(k=1;k<=l;k++) scale += fabs(a[i][k]);
      /* малая величина шкалы -> пропустить преобразование */
      if(scale==0.) e[i]=a[i][l];
      else {
       /* отмасштабировать строку и вычислить s2 в h */
        for(k=1;k<=l;k++) {
          a[i][k]/=scale; h += a[i][k]*a[i][k];
        }
       /* вычислить вектор u */
        f=a[i][l];
        g=(f>=0.?-sqrt(h):sqrt(h));
        e[i]=scale*g; h -= f*g;
       /* записать u на место i-го ряда a */
        a[i][l]=f-g;
       /* вычисление u/h, Au, p, K */
        f=0.;
        for(j=1;j<=l;j++) {
       /* следующая инструкция не нужна, если не требуются вектора,
             она содержит загрузку u/h в столбец a */
          a[j][i]=a[i][j]/h;
       /* сформировать элемент Au (в g) */
          g=0.;
          for(k=1;k<=j;k++) g += a[j][k]*a[i][k];
          for(k=j+1;k<=l;k++) g += a[k][j]*a[i][k];
       /* загрузить элемент p во временно неиспользуемую область e */
          e[j]=g/h;
       /* подготовка к формированию K */
          f += e[j]*a[i][j];
        }
       /* Сформировать K */
        hh=f/(h+h);
        for(j=1;j<=l;j++) {
       /* Сформировать q и поместить на место p (в e) */
          f=a[i][j]; e[j]=g=e[j]-hh*f;
       /* Трансформировать матрицу a */
          for(k=1;k<=j;k++) a[j][k] -= (f*e[k]+g*a[i][k]);
        }
      }
    }
    else e[i]=a[i][l];
    d[i]=h;
  }
  /* если не нужны собственные вектора, опустите следующую инструкцию */
  d[1]=0.;
  /* эту опускать не надо */
  e[1]=0.;
  /* Все содержание цикла, кроме одной инструкции, можно опустить, если не
     требуются собственные вектора */
  for(i=1;i<=n;i++) {
    l=i-1;
    /* этот блок будет пропущен при i=1 */
    if(d[i]!=0.) {
      for(j=1;j<=l;j++) {
        g=0.;
       /* формируем PQ, используя u и u/H */
        for(k=1;k<=l;k++) g += a[i][k]*a[k][j];
        for(k=1;k<=l;k++) a[k][j] -= g*a[k][i];
      }
    }
    /* эта инструкция остается */
    d[i]=a[i][i];
    /* ряд и колонка матрицы a преобразуются к единичной, для след. итерации */
    a[i][i]=0.;
    for(j=1;j<=l;j++) a[j][i]=a[i][j]=0.;
  }
}