当前位置：首页 > news >正文

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

news 来源：原创 2024/9/20 5:08:52

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_decall: $(EXE)%: %.cpp utils.cppg++ $^ -o $@ -lm -DBUILD_MAINhouseholder_qr_dec: householder_qr_dec.cpp householder_vector.cpp utils.cppg++ -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 )^tif(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);