import numpy as np

n = 100
a = np.random.rand(n)
b = a + 1e-9 * np.random.randn(n)

# Сначала нормируем вектор a
a = a / np.linalg.norm(a)

# Убираем из b проекцию на a (попытка сделать c ортогональным к a)
c = b - np.dot(b, a) * a
c = c / np.linalg.norm(c)

# Скалярное произведение c и a пока не идеально мало:
print(np.dot(c, a))

# Делаем повторную ортогонализацию:
# Снова вычитаем проекцию c на a
c = c - np.dot(c, a) * a
c = c / np.linalg.norm(c)

# Теперь скалярное произведение гораздо меньше!
print(np.dot(c, a))

-2.7671064037337878e-08
-6.591949208711867e-17

Упражнение

В этом упражнении мы реализуем классический и модифицированный процессы Грама-Шмидта для построения QR-разложения матрицы.

Вход: матрица A \in \mathbb{R}^{n \times n} Выход: матрицы Q и R такие, что A = QR, при этом Q ортогональная, а R верхнетреугольная. Мы сфокусируемся на матрице Q, потому что имея её, мы можем легко получить R по формуле R = Q^T A.

• Численная неустойчивость. Классический процесс Грама-Шмидта (CGS) может приводить к большим погрешностям при вычислении ортонормированной системы векторов, особенно при плохой обусловленности исходной матрицы. В машинной арифметике возникает потеря ортогональности.

• Модифицированный процесс (MGS). Идея заключается в том, чтобы последовательно вычитать проекции на уже построенные ортонормированные вектора, шаг за шагом. Формально, если q_1, \dots, q_{k-1} уже построены, то для k-го вектора: q_k := a_k, \quad q_k := q_k - (q_k,\,q_1)\,q_1, \quad q_k := q_k - (q_k,\,q_2)\,q_2, \quad \dots В точной арифметике это эквивалентно классическому процессу, но в машинной арифметике MGS, как правило, даёт более устойчивое решение.

• Сложность. Оба алгоритма имеют асимптотическую сложность порядка O(n^2 m), где m \times n — размерность исходной матрицы. Для квадратной матрицы n \times n сложность O(n^3).

• QR-разложение. По готовой ортонормированной матрице Q можно найти верхнетреугольную R, например, умножив Q^T A. Тогда A = Q R.

import numpy as np

def hilbert_matrix(n):
    H = np.zeros((n, n), dtype=float)
    for i in range(n):
        for j in range(n):
            H[i, j] = 1.0 / (i + j + 1)
    return H

def generate_matrix(n, type='random'):
    if type == 'random':
        A = np.random.randn(n, n)
    elif type == 'hilbert':
        A = hilbert_matrix(n)
    return A

def classical_gs(A):
    m, n = A.shape
    Q = np.zeros((m, n), dtype=float)
    for k in range(n):
        # Берём исходный столбец
        v = A[:, k].copy()
        
        # Считаем все скалярные произведения \alpha_j = (q_j, v) заранее
        # (здесь q_j уже ортонормированы из предыдущих шагов)
        alphas = np.array([np.dot(Q[:, j], v) for j in range(k)])
        
        # Вычитаем сумму \sum_j \alpha_j * q_j
        for j in range(k):
            v -= alphas[j] * Q[:, j]
        
        # Нормируем
        norm_v = np.linalg.norm(v)
        Q[:, k] = v / norm_v
    return Q

def modified_gs(A):
    m, n = A.shape
    Q = np.zeros((m, n), dtype=float)
    for k in range(n):
        v = A[:, k].copy()
        # В МОДИФИЦИРОВАННОМ алгоритме на каждом шаге j
        # берём текущее "свежее" v
        for j in range(k):
            alpha = np.dot(Q[:, j], v)
            v -= alpha * Q[:, j]
        Q[:, k] = v / np.linalg.norm(v)
    return Q


# Пример: проверяем ортогонализацию на случайной матрице
n = 10
A = generate_matrix(n, type='hilbert')

Q_cgs = classical_gs(A)
Q_mgs = modified_gs(A)

# Оцениваем "ошибку ортогональности" как ||Q^T Q - I|| по Фробениусу
err_cgs = np.linalg.norm(Q_cgs.T @ Q_cgs - np.eye(n), ord='fro')
err_mgs = np.linalg.norm(Q_mgs.T @ Q_mgs - np.eye(n), ord='fro')

# Фробениусова норма единичной матрицы I (это sqrt(n))
norm_I_fro = np.linalg.norm(np.eye(n), ord='fro')

rel_cgs_fro = err_cgs / norm_I_fro
rel_mgs_fro = err_mgs / norm_I_fro

print("Абсолютная ошибка (CGS) =", err_cgs)
print("Абсолютная ошибка (MGS) =", err_mgs)
print("Относительная ошибка (CGS) =", rel_cgs_fro)
print("Относительная ошибка (MGS) =", rel_mgs_fro)

Абсолютная ошибка (CGS) = 3.4653383412436143
Абсолютная ошибка (MGS) = 8.763548618214397e-05
Относительная ошибка (CGS) = 1.095836202143963
Относительная ошибка (MGS) = 2.7712774019178854e-05