基于 Householder 变换的 qr 分解算法与源码实现

news2024/9/20 0:55:29

1，算法描述

1.1 算法1 反射向量

计算 Householder 向量

给定 $\mathbf{x} \in \mathbb{R}^m$ 算法计算满足 v(1) = 1.0 的 $\mathbf{v}\in \mathbb{R}^m$ 和 $\beta \in \mathbb{R}$ , 使得 $\mathbf{P} = \mathbf{I}_m - \beta \mathbf{v}\mathbf{v}^ T$ 是正交矩阵且 $\mathbf{Px} = \parallel \mathbf{x} \parallel_2\mathbf{e}_1$ , 即，将m维向量 $\mathbf{x}$ 通过反射变换 $\mathbf{P}$ 反射至 $\mathbf{e}_1$ 轴上去。

$\mathbf{function} [v, \beta] = \mathbf{house(x)}$

$m = \mathbf{length(x)}$

$\sigma = \mathbf{x(2:m)^T x(2:m)}$

$\mathbf{v}=\left[\begin{array}{c} 1\\ \mathbf{x(2:m)}\end{array} \right ]$

$\mathbf{if( } \sigma==0 \&\& \mathbf{x(1) >= 0} \mathbf{)}$

$\beta = 0$

$\mathbf{elseif(} \sigma==0 \&\& \mathbf{x(1)} <0 \mathbf{)}$

$\beta = 2.0$

$\mathbf{else}$

$u = \sqrt{\mathbf{x(1)}+\sigma}$

$\mathbf{if(x(1)}<=0\mathbf{)}$

$\mathbf{v(1) = x(1)} - u$

$\mathbf{else}$

$\mathbf{v(1)} = -\sigma/(\mathbf{x(1)}+u)$

$\mathbf{end}$

$\beta = 2*\mathbf{v(1)}^2/(\sigma + \mathbf{v(1)}^2)$

$\mathbf{v} = \mathbf{v}/\mathbf{v(1)}$

$\mathbf{end}$

1.2 算法2 QR 分解

Householder QR 分解

未完待补。。。。

2，源码实现

由如下几个文件组成：

Makefile
householder_qr_dec.cpp
householder_vector.cpp
householder_vector.h
utils.cpp
utils.h

Makefile


EXE := householder_qr_dec

all: $(EXE)

%: %.cpp utils.cpp
	g++ $^ -o $@ -lm -DBUILD_MAIN

householder_qr_dec: householder_qr_dec.cpp householder_vector.cpp utils.cpp
	g++ -g $^ -o $@ -lm -DBUILD_MAIN

.PHONY: clean
clean:
	-rm -rf $(EXE)

householder_qr_dec.cpp

#include <stdlib.h>
#include <string.h>
#include "utils.h"
#include "householder_vector.h"

void vector_mul_matrix(int m, int n, float* vt, float* A, int lda, float* C)// C(1, n) = vt(1, m) * A(m, n)
{
    for(int j=0; j<n; j++){
        float sigma = 0.0f;
        for(int k=0; k<m; k++){
            sigma += vt[k] * A[k + j*lda];
        }
        C[j] = sigma;
    }
}

// A(m, n) = A - b * v(m) * C(n)
void rank_one_update(int m, int n, float beta, float* v, float* C, float* A, int lda)
{
    for(int i=0; i<m; i++){
        for(int j=0; j<n; j++){
            A[i + j*lda] -= beta * v[i] * C[j];
        }
    }
}

void store_householder_vector(float* A, float* v, int len)
{
    for(int i=0; i<len; i++){
        A[i] = v[i];
    }
}

// (Householder QR)
void householder_qr(int m, int n, float* A, int lda, float* tau)
{
    float* v = nullptr;
    float beta = 0;
    float* C = nullptr;

    C = (float*)malloc(m*sizeof(float));
    v = (float*)malloc(m*sizeof(float));

    for(int j=0; j<n; j++){
        //step1, [v, beta] = house(A(j:m, j))
        //void house(int m, float* x, float* v, float& beta);
        house(m-j, A+j+j*lda, v, beta);
        tau[j] = beta;
        //step2, A(j:m, j:n) = (I - beta*v*v^T)A(j:m, j:n)//m 可能挺大，n<=32 or 64;
        // A = A - b*v*(v^t * A) = A - b*v*C;  (j:m, j:n),   v(j:m),   (v^t * A)(j:n)

          //step2.1 C(j:n) = (v^t)(1, j:m) * A(j:m, j:n)
        vector_mul_matrix(m-j, n-j, v, A+j+j*lda, lda, C);

          //step2.2 A(j:m, j:n) = A -  b*v*C;
        rank_one_update(m-j, n-j, beta, v, C, A+j+j*lda, lda);
        //step3, if(j<m) A(j+1 : m, j) = v(2 : m-j+1)
        store_householder_vector(A+(j+1 + j*lda), v+1, m-j-1);
    }
}




int main()
{
    int m = 4;
    int n = 4;
    int lda = m;
    int ldb = lda;
    float* A = nullptr;
    float* B = nullptr;

    A = (float*)malloc(lda*n*sizeof(float));
    B = (float*)malloc(ldb*n*sizeof(float));

    init_matrix(m, n, A, lda, 2024); printf("A =[ ...\n");   print_matrix(m, n, A, lda);
    memcpy(B, A, lda*n*sizeof(float));

    float* tau = nullptr;
    tau = (float*)malloc(n*sizeof(float));

    householder_qr(m, n, A, lda, tau);
    printf("R+tau =\n");
    print_matrix(m, n, A, lda);

    printf("\ntau = \n");print_vector(n, tau);

    return 0;
}

householder_vector.cpp

#include <stdlib.h>


#include "utils.h"

//#define BUILD_MAIN
void v_is_1_x_2_m(int M, float* x, float* v)
{
    v[0] = 1.0f;
    for(int i=1; i<M; i++){
        v[i] = x[i];
    }
}

void dot_vector(int n, float* x, float* y, float& sigma)
{
    sigma = 0.0f;

    for(int i=0; i<n; i++)
    {
        sigma += x[i]*y[i];
    }
}

void vector_div_scalar(int M, float* v, float alpha)
{
    for(int i=0; i<M; i++){
        v[i] /= alpha;
    }
}

void house(int m, float* x, float* v, float& beta)
{
    float sigma;
    dot_vector(m-1, x+1, x+1, sigma);//    printf("in house() sigma = %7.3f\n", sigma);
    v_is_1_x_2_m(m, x, v);// v= ( 1 x(2 : m)^t )^t

    if(sigma==0.0f && x[0]>=0.0f)
    {
        beta = 0.0f;
    }
    else if(sigma==0.0f && x[0]<0.0f)
    {
        beta = 2.0f;
    }
    else
    {
        float miu;

        miu = sqrt(x[0]*x[0] + sigma);
        if(x[0]<= 0.0f){
            v[0] = x[0] - miu;
        }
        else{
            v[0] = -sigma/(x[0]+miu);
        }

        beta = 2.0f*v[0]*v[0]/(sigma + v[0]*v[0]);
        vector_div_scalar(m, v, v[0]);
    }
}

householder_vector.h

#pragma once 

void house(int m, float* x, float* v, float& beta);

utils.cpp

#include "utils.h"

void init_matrix(int M, int N, float* A, int lda, int seed)
{
    srand(seed);

    for(int i=0; i<M; i++){
        for(int j=0; j<N; j++){

            int r;
            r = rand();
            A[i + j*lda] = (((r>(RAND_MAX/3)) ? 3.0f : -3.0f)*rand())/RAND_MAX;
        }
    }
}

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

void print_vector(int N, float* A)
{
    print_matrix(1, N, A, 1);
}

utils.h

#pragma once
#include <stdio.h>
#include <math.h>

void init_matrix(int M, int N, float* A, int lda, int seed);
void print_matrix(int M, int N, float* A, int lda);
void print_vector(int M, float* A);