基于Givens矩阵的QR矩阵分解

QR分解是一种将矩阵分解为正交矩阵和上三角矩阵的方法。在QR分解中，正交矩阵Q的转置是它的逆矩阵，因此QR分解可以用于求解线性方程组、最小二乘问题等。

二阶Givens矩阵

一般地，二阶Givens矩阵记为下列形式：

其中

下面开始介绍基于Givens矩阵的QR分解算法。Givens矩阵是一种旋转矩阵，可以将一个向量旋转到另一个向量的方向。在QR分解中，我们使用Givens矩阵将矩阵的列向量逐个旋转，使得矩阵变为上三角矩阵。

QR分解的详细步骤如下：

1. 对矩阵A的第一列进行Givens变换，使得A的第一列的下面的元素都变为0。这样得到一个新的矩阵A1和一个Givens矩阵G1。

2. 对矩阵A1的第二列进行Givens变换，使得A1的第二列的下面的元素都变为0。这样得到一个新的矩阵A2和一个Givens矩阵G2。

3. 重复步骤2，直到所有的列都被处理完毕。这样得到一个上三角矩阵R和一个正交矩阵Q，满足A=QR。

大致如下：

下面是基于Givens矩阵的QR分解的C和Python代码：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define N 3 // 指定矩阵阶数

void print_matrix(double **M) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            printf("%f ", M[i][j]);
        }
        printf("\n");
    }
}

// 矩阵乘法 
double **matrix_multiply(double **A, double **B) {
    int i, j, k;
    double **C = (double **)malloc(sizeof(double *) * N);
    for (i = 0; i < N; i++) {
        C[i] = (double *)malloc(sizeof(double) * N);
    }
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            C[i][j] = 0;
            for (k = 0; k < N; k++) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }

    return C;
}

// 单位矩阵 
double **identity_matrix() {
    int i, j;
    double **I = (double **)malloc(sizeof(double *) * N);
    for (i = 0; i < N; i++) {
        I[i] = (double *)malloc(sizeof(double) * N);
        for (j = 0; j < N; j++) {
            if (i == j) {
                I[i][j] = 1;
            } else {
                I[i][j] = 0;
            }
        }
    }
    
    return I;
}

// 矩阵的转置
double **transpose_matrix(double **A) {
    int i, j;
    double **AT = (double **)malloc(sizeof(double *) * N);
    for (i = 0; i < N; i++) {
        AT[i] = (double *)malloc(sizeof(double) * N);
        for (j = 0; j < N; j++) {
            AT[i][j] = A[j][i];
        }
    }
    
    return AT;
}

// Givens变换c和s 
double *givens_rotation(double a, double b) {
    double c, s;
    double *cs = (double *)malloc(sizeof(double) * 2);
    if (b == 0) {
        c = 1;
        s = 0;
    } else if (fabs(b) > fabs(a)) {
        double r = a / b;
        s = 1 / sqrt(1 + r * r);
        c = s * r;
    } else {
        double r = b / a;
        c = 1 / sqrt(1 + r * r);
        s = c * r;
    }
    cs[0] = c;
    cs[1] = s;

    return cs;
}

double **qr_givens(double A[][N], double **Q) {
    double **Q_ = identity_matrix();
    Q = Q_;
    double **R = (double **)malloc(sizeof(double *) * N);
    for (int i = 0; i < N; i++) {
        R[i] = (double *)malloc(sizeof(double) * N);
    }
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            R[i][j] = A[i][j];
        }
    }

    for (int j = 0; j < N; j++) {
        for (int i = N - 1; i > j; i--) {
            double **G = identity_matrix();
            double *cs = givens_rotation(R[i - 1][j], R[i][j]);
            double c = cs[0];
            double s = cs[1];
            G[i - 1][i - 1] = c;
            G[i - 1][i] = s;
            G[i][i - 1] = -s;
            G[i][i] = c;
            
            double **R_ = matrix_multiply(G, R);
            R = R_;
            double **GT = transpose_matrix(G);
            double **Q_ = matrix_multiply(Q, GT);
            Q = Q_;
        }
    }

    return R;
}

int main() {
    double A[N][N] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
    double **Q = (double **)malloc(sizeof(double *) * N);
    for (int i = 0; i < N; i++) {
        Q[i] = (double *)malloc(sizeof(double) * N);
    }

    double **R = qr_givens(A, Q);
    print_matrix(R);
    
    return 0;
}

import numpy as np


def givens_rotation(a, b):  # 定义一个Givens旋转函数，a和b是两个数
    if b == 0:  # 如果b为0
        return 1, 0  # 返回1和0
    elif abs(b) > abs(a):  # 如果b的绝对值大于a的绝对值
        r = a / b  # r等于a除以b
        s = 1 / np.sqrt(1 + r ** 2)  # s等于1除以根号下1加r的平方
        c = s * r  # c等于s乘以r
    else:  # 否则
        r = b / a  # r等于b除以a
        c = 1 / np.sqrt(1 + r ** 2)  # c等于1除以根号下1加r的平方
        s = c * r  # s等于c乘以r
    return c, s  # 返回c和s


def qr_givens(A):  # 定义一个QR分解函数，A是一个矩阵
    m, n = A.shape  # m和n分别等于A的行数和列数
    Q = np.identity(m)  # Q等于m阶单位矩阵
    R = A.copy()  # R等于A的副本
    for j in range(n):  # 对于j在0到n-1的范围内(从第一列到第n列)
        for i in range(m - 1, j, -1):  # 对于m-1到j+1的范围内(从最后一行到第j行)
            G = np.identity(m)  # G等于m阶单位矩阵
            c, s = givens_rotation(R[i - 1, j], R[i, j])  # c和s等于Givens旋转函数的返回值
            G[i - 1:i + 1, i - 1:i + 1] = [[c, s], [-s, c]]  # G的第i-1到i行，第i-1到i列等于[[c, s], [-s, c]]
            R = np.dot(G, R)  # R等于G和R的矩阵乘积
            Q = np.dot(Q, G.T)  # Q等于Q和G的转置矩阵的矩阵乘积
    return Q, R  # 返回Q和R


A = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])  # A等于一个3行3列的矩阵
Q, R = qr_givens(A)  # Q和R等于QR分解函数的返回值
print("Q:\n", Q)  # 输出Q