数值分析作业 - 矩阵特征值的数值解法

数值分析

Problem 1

幂迭代: 找出绝对值最大的特征值, 为
逆向幂迭代: 找出离最近的特征值, 为
幂迭代的线性收敛率满足 ; 逆向幂迭代先平移为再取倒数此时, 则有 , 逆向幂迭代收敛更快.

Problem 2

约化第一列

约化第二列

则有

Problem 3

使用 PyTorch 实现

幂迭代

from typing import List, Tuple
import torch
import math

def power_method(A: torch.Tensor, x0: torch.Tensor, limit: int, tol=1e-6) -> List[Tuple[float, torch.Tensor]]:
    """
    Power method for computing the largest eigenvalue of a matrix.
    :param A: A square matrix
    :param x0: An initial vector
    :param limit: The number of iterations
    :param tol: The tolerance for convergence
    :return: A list of tuples (eigenvalue, eigenvector)
    """
    x = x0
    for k in range(limit):
        i = x.abs().argmax()
        x = x / x[i]
        y = A @ x
        t = y[i].item()
        e1 = (y / t - x).norm()
        if e1 < tol:
            return [(t, y)]
        z = A @ y
        t = z[i].item()
        e2 = (z / t - x).norm()
        if e2 < tol and e1 > tol * 100.:
            lam = math.sqrt(t)
            return [(lam, z + lam * y), (-lam, z - lam * y)]
        x = z
    return []

逆幂迭代

def inverse_power_method(A: torch.Tensor, x0: torch.Tensor, limit: int, tol=1e-6) -> List[Tuple[float, torch.Tensor]]:
    """
    Inverse power method for computing the smallest eigenvalue of a matrix.
    :param A: A square matrix
    :param x0: An initial vector
    :param limit: The number of iterations
    :param tol: The tolerance for convergence
    :return: A list of tuples (eigenvalue, eigenvector)
    """
    x = x0
    LU, pivots = torch.linalg.lu_factor(A)
    for k in range(limit):
        i = x.abs().argmax()
        x = x / x[i]
        y = torch.linalg.lu_solve(LU, pivots, x.view(-1, 1)).view(-1)
        t = y[i].item()
        e1 = (y / t - x).norm()
        if e1 < tol:
            return [(1 / t, x)]
        z = torch.linalg.lu_solve(LU, pivots, y.view(-1, 1)).view(-1)
        t = z[i].item()
        e2 = (z / t - x).norm()
        if e2 < tol and e1 > tol * 100.:
            lam = 1 / math.sqrt(t)
            return [(lam, x + lam * y), (-lam, x - lam * y)]
        x = z
    return []

运行算法求解矩阵

A = torch.tensor([[8, -8, -4], [12, -15, -7], [-18, 26, 12]], dtype=torch.float64)
x0 = torch.tensor([1, 1, 1], dtype=torch.float64)
torch.set_printoptions(precision=8, sci_mode=False)
print("Using the power method:")
print(power_method(A, x0, 12))
print("Using the inverse power method:")
print(inverse_power_method(A, x0, 12))

输出

Using the power method:
[(3.9999996821089994, tensor([-3.99999873, -3.99999905,  3.99999968], dtype=torch.float64))]
Using the inverse power method:
[(-0.9999998211860279, tensor([    -0.00000024,     -0.50000018,      1.00000000],
       dtype=torch.float64))]

接近理论结果.

Problem 4

QR 分解

def qr_decomposition_full(A: torch.Tensor, eps=1e-8) -> Tuple[torch.Tensor, torch.Tensor]:
    """
    Full QR decomposition of a matrix.
    :param A: A matrix
    :param eps: The tolerance for convergence
    :return: A tuple (Q, R) where Q is an orthogonal matrix and R is an upper triangular matrix
    """
    m, n = A.shape
    Q = torch.eye(m, m, dtype=A.dtype)
    R = A.clone()
    for k in range(n - 1):
        t = R[k:, k].norm()
        if t < eps:
            continue
        if R[k, k] > 0:
            t = -t
        v = -R[k:, k].clone()
        v[0] += t
        R[k, k] = t
        R[k + 1:, k] = 0
        t0 = 2 / (v @ v)
        for j in range(k + 1, n):
            t = t0 * (v @ R[k:, j])
            R[k:, j] -= t * v
        for i in range(m):
            t = t0 * (Q[i, k:] @ v)
            Q[i, k:] -= t * v
    return Q, R

QR 算法

def qr_method(A: torch.Tensor, limit: int) -> Tuple[torch.Tensor, torch.Tensor]:
    """
    QR method for computing the eigenvalues of a matrix.
    :param A: A square matrix
    :param limit: The number of iterations
    :param tol: The tolerance for convergence
    :return: A tuple (eigenvalues, Q) where eigenvalues is a vector of eigenvalues and Q is an orthogonal matrix
    """
    n = A.shape[0]
    Qbar = torch.eye(n, n, dtype=A.dtype)
    for k in range(limit):
        Q, R = qr_decomposition_full(A)
        A = R @ Q
        Qbar = Qbar @ Q
    return A.diag(), Qbar

运行求解

1 2	print("Using the QR method:") print(qr_method(A, 25))

输出:

Using the QR method:
(tensor([ 4.00000001,  1.99999996, -0.99999997], dtype=torch.float64), 
tensor([[-0.57735026, -0.61721341,  0.53452248],
        [-0.57735027, -0.15430336, -0.80178373],
        [ 0.57735028, -0.77151674, -0.26726124]], dtype=torch.float64))

注意到只有特征值的特征向量收敛到了正确结果

本文采用署名-相同方式共享 4.0 国际许可协议，转载请注明出处。