Householder transformation + QL to calculate eigenValue and eigenVectors of Hertian Matrix, cpu code

article/2025/10/3 8:03:18

化Hertinan matrix eigen problem into a real symmetric matrix eigen problem:

原理：

与matlab的结果进行对比：

A=[ ...
( 3.0 + 0.0*j) (-2.0   -2.0*j) (-0.9    -0.9*j) (-0.5    -0.5*j); ...
(-2.0 + 2.0*j) ( 4.0 +  0.0*j) ( 1.0 +   1.0*j) (-0.7    -0.5*j); ...
(-0.9 + 0.9*j) ( 1.0   -1.0*j) (-1.0 +   0.0*j) ( 0.1 +   0.1*j); ...
(-0.5 + 0.5*j) (-0.7 +  0.5*j) ( 0.1    -0.1*j) ( 1.0 +   0.0*j); ...
][V D]=eig(A)

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>using namespace std;
struct float2 {float x;float y;
};typedef struct float2  float2;void print_matrix(float* A, int M, int N, int lda);
void print_matrix_complex(float2* A, int M, int N, int lda);
void print_vector(float* A, int N);
/
void real_neimag_real_imag_fill(float2* A_h, int N, int lda, float* S, int n, int lds);
float max(float x, float y);
void swap(float& x, float& y);
void tred2(float* V, float* d, float* e, int n, int ldv);
void tql2(float* V, float* d, float* e, int n, int ldv);
void eigen_vector_decompose_hermitian(float2* A, int N, int lda, float2* V_cevd, float* d_cevd);
/#define NA 4/*
* complex A=(4x4)  math topology:( 3.0,  0.0),  (-2.0,  -2.0),  (-0.9,  -0.9),  (-0.5,  -0.5),(-2.0,  2.0),  ( 4.0,   0.0),  ( 1.0,   1.0),  (-0.7,  -0.5),(-0.9,  0.9),  ( 1.0,  -1.0),  (-1.0,   0.0),  ( 0.1,   0.1),(-0.5,  0.5),  (-0.7,   0.5),  ( 0.1,  -0.1),  ( 1.0,   0.0)float A_h[NA * (2*NA)] = {//complex columnMajor: mem topology3.0, 0.0,  -2.0, 2.0,  -0.9, 0.9,  -0.5,0.5,-2.0,-2.0,   4.0, 0.0,   1.0,-1.0,  -0.7, 0.5,-0.9,-0.9,   1.0, 1.0,  -1.0, 0.0,   0.1,-0.1,-0.5,-0.5,  -0.7,-0.5,   0.1, 0.1,   1.0, 0.0*/int main() {float A_h[NA * (2 * NA)] = {//complex columnMajor:3.0, 0.0,  -2.0, 2.0,  -0.9, 0.9,  -0.5, 0.5,-2.0,-2.0,   4.0, 0.0,   1.0,-1.0,  -0.7, 0.5,-0.9,-0.9,   1.0, 1.0,  -1.0, 0.0,   0.1,-0.1,-0.5,-0.5,  -0.7,-0.5,   0.1, 0.1,   1.0, 0.0};///Parameters:int N = NA;int lda = N;float2* A = nullptr;
//___________________________________________________________________A = (float2*)malloc(lda * N * sizeof(float2));memcpy(A, A_h, lda * N * sizeof(float2));printf("A(ma)=\n");print_matrix_complex(A, N, N, N);
//______________________________________________________________//    void eigen_vector_decompose(float2 * A, int N, int lda, float2 * V_cevd, float* d_cevd);float* d_cevd = nullptr;float2* V_cevd = nullptr;int ldv_cevd = N;d_cevd = (float*)malloc(N * sizeof(float));V_cevd = (float2*)malloc(ldv_cevd * N * sizeof(float2));eigen_vector_decompose_hermitian(A, N, lda, V_cevd, d_cevd);#if 0int n = 2 * N;int lds = n;float* S = nullptr;S = (float*)malloc(lds * n * sizeof(float));real_neimag_real_imag_fill((float2*)A, N, lda, S, n, lds);printf("S(ma)=\n");print_matrix(S, n, n, lds);
//______________________________________________________________int ldv_evd = n;float* V_evd = nullptr;V_evd = (float*)malloc(lds * n * sizeof(float));memcpy(V_evd, S, lds * n * sizeof(float));printf("V_evd(ma)=\n");print_matrix(V_evd, n, n, ldv_evd);float* e = nullptr;float* d_evd = nullptr;e = (float*)malloc(n * sizeof(float));d_evd = (float*)malloc(n * sizeof(float));//______________________________________________________________//void tred2(float* V, float* e, float* d, int n, int ldv);tred2(V_evd, d_evd, e, n, ldv_evd);printf("V_evd-tred2=\n");print_matrix(V_evd, n, n, ldv_evd);printf("\nd-tred2=\n");print_vector(d_evd, n);printf("\ne-tred2=\n");print_vector(e, n);//void tql2(float* V, float* d, float* e, int n, int ldv);tql2(V_evd, d_evd, e, n, ldv_evd);printf("\nV_evd-tql2=\n");print_matrix(V_evd, n, n, ldv_evd);printf("\nd-tql2=\n");print_vector(d_evd, n);printf("\ne-tql2=\n");print_vector(e, n);
//_________________________________________________________________________________/**float* d_cevd = nullptr;float2* V_cevd = nullptr;int ldv_cevd = N;d_cevd = (float*)malloc(N * sizeof(float));V_cevd = (float2*)malloc(ldv_cevd * N * sizeof(float2));*/for(int i=0; i<N; i++)d_cevd[i] = d_evd[2 * i];printf("\nd_cevd=\n");print_vector(d_cevd, N);for (int j = 0; j < N; ++j){int j2 = 2 * j;for (int i = 0; i < N; ++i) {V_cevd[i + j * ldv_cevd].x = V_evd[i + j2 * ldv_evd];//complex<Type>(RV[i][j2], RV[i + N][j2]);V_cevd[i + j * ldv_cevd].y = V_evd[(i + N) + j2 * ldv_evd];}}printf("\nV_cevd=\n");print_matrix_complex(V_cevd, N, N, ldv_cevd);#endifreturn 0;
}void real_neimag_real_imag_fill( float2 * A_h, int N, int lda, float* S, int n, int lds)//n=2*N
{float2 *A;A = A_h;for (int i = 0; i < N; i++) {for (int j = 0; j < N; j++) {S[i +j*lds] = A[i+ j*lda].x;S[i +(j + N)*lds] = -A[i +j*lda].y;S[(i + N)+ j*lds] = -S[i +(j + N)*lds];S[(i + N)+ (j + N)*lds] = S[i+ j*lds];}}
}//typedef float float;void tred2(float *V, float* d, float *e, int n, int ldv)
{for (int j = 0; j < n; ++j)//LL::d[j] = V[n - 1][j];d[j] = V[(n - 1)+j*ldv];// Householder reduction to tridiagonal formfor (int i = n - 1; i > 0; --i){// scale to avoid under/overflowfloat scale = 0;float h = 0;for (int k = 0; k < i; ++k) {//   scale += fabs(d[k]);printf(" d[k]=%f", d[k]);scale += fabs(d[k]);printf(" scale=%f", scale);}printf("\n______________\n");if (scale == 0){e[i] = d[i - 1];for (int j = 0; j < i; ++j){d[j] = V[(i - 1)+j*ldv];V[i+j*ldv] = 0;V[j+i*ldv] = 0;}}else{// generate Householder vectorfor (int k = 0; k < i; ++k){d[k] /= scale;h += d[k] * d[k];}float f = d[i - 1];float g = sqrt(h);if (f > 0)g = -g;e[i] = scale * g;h = h - f * g;d[i - 1] = f - g;for (int j = 0; j < i; ++j)e[j] = 0;// Apply similarity transformation to remaining columns.for (int j = 0; j < i; ++j){f = d[j];V[j+i*ldv] = f;g = e[j] + V[j+j*ldv] * f;for (int k = j + 1; k <= i - 1; ++k){g += V[k+j*ldv] * d[k];e[k] += V[k+j*ldv] * f;}e[j] = g;}f = 0;for (int j = 0; j < i; ++j){e[j] /= h;f += e[j] * d[j];}float hh = f / (h + h);for (int j = 0; j < i; ++j)e[j] -= hh * d[j];for (int j = 0; j < i; ++j){f = d[j];g = e[j];for (int k = j; k <= i - 1; ++k)V[k+j*ldv] -= (f * e[k] + g * d[k]);d[j] = V[(i - 1)+j*ldv];V[i+j*ldv] = 0;}}d[i] = h;}// accumulate transformationsfor (int i = 0; i < n - 1; i++){V[n - 1+i*ldv] = V[i+i * ldv];V[i+i * ldv] = 1;float h = d[i + 1];if (h != 0){for (int k = 0; k <= i; ++k)d[k] = V[k+(i + 1) *ldv] / h;for (int j = 0; j <= i; ++j){float g = 0;for (int k = 0; k <= i; ++k)g += V[k+(i + 1) * ldv] * V[k+j * ldv];for (int k = 0; k <= i; ++k)V[k+j * ldv] -= g * d[k];}}for (int k = 0; k <= i; ++k)V[k+(i + 1) * ldv] = 0;}for (int j = 0; j < n; ++j){d[j] = V[(n - 1)+j * ldv];V[(n - 1)+j * ldv] = 0;}V[(n - 1)+(n - 1) * ldv] = 1;e[0] = 0;
}float max(float x, float y) {return (x > y) ? x : y;
}void swap(float& x, float& y) {float tmp;tmp = x;x = y;y = tmp;
}void tql2(float* V, float* d, float* e, int n, int ldv)
{for (int i = 1; i < n; ++i)e[i - 1] = e[i];e[n - 1] = 0;float f = 0;float tst1 = 0;float eps = pow(2.0, -52.0);for (int l = 0; l < n; ++l){// find small subdiagonal elementtst1 = max(tst1, fabs(d[l]) + fabs(e[l]));int m = l;// original while-loop from Java codewhile (m < n){if (fabs(e[m]) <= eps * tst1)break;m++;}// if m == l, d[l] is an eigenvalue, otherwise, iterateif (m > l){int iter = 0;do{iter = iter + 1;// compute implicit shiftfloat g = d[l];float p = (d[l + 1] - g) / (2.0 * e[l]);float r = hypot(p, 1.0);if (p < 0)r = -r;d[l] = e[l] / (p + r);d[l + 1] = e[l] * (p + r);float dl1 = d[l + 1];float h = g - d[l];for (int i = l + 2; i < n; ++i)d[i] -= h;f += h;// implicit QL transformation.p = d[m];float c = 1;float c2 = c;float c3 = c;float el1 = e[l + 1];float s = 0;float s2 = 0;for (int i = m - 1; i >= l; --i){c3 = c2;c2 = c;s2 = s;g = c * e[i];h = c * p;r = hypot(p, e[i]);e[i + 1] = s * r;s = e[i] / r;c = p / r;p = c * d[i] - s * g;d[i + 1] = h + s * (c * g + s * d[i]);// accumulate transformation.for (int k = 0; k < n; ++k){h = V[k+(i + 1)*ldv];V[k+(i + 1) * ldv] = s * V[k+i * ldv] + c * h;V[k+i * ldv] = c * V[k+i * ldv] - s * h;}}p = -s * s2 * c3 * el1 * e[l] / dl1;e[l] = s * p;d[l] = c * p;} while (fabs(e[l]) > eps * tst1);}d[l] += f;e[l] = 0;}// Sort eigenvalues and corresponding vectors.for (int i = 0; i < n - 1; ++i){int k = i;float p = d[i];for (int j = i + 1; j < n; ++j)if (d[j] < p){k = j;p = d[j];}if (k != i){d[k] = d[i];d[i] = p;for (int j = 0; j < n; ++j)swap(V[j+i * ldv], V[j+k * ldv]);}}
}void eigen_vector_decompose_hermitian(float2* A, int N, int lda, float2* V_cevd, float* d_cevd)
{int n = 2 * N;int lds = n;float* S = nullptr;S = (float*)malloc(lds * n * sizeof(float));real_neimag_real_imag_fill((float2*)A, N, lda, S, n, lds);printf("S(ma)=\n");print_matrix(S, n, n, lds);//______________________________________________________________int ldv_evd = n;float* V_evd = nullptr;V_evd = (float*)malloc(lds * n * sizeof(float));memcpy(V_evd, S, lds * n * sizeof(float));printf("V_evd(ma)=\n");print_matrix(V_evd, n, n, ldv_evd);float* e = nullptr;float* d_evd = nullptr;e = (float*)malloc(n * sizeof(float));d_evd = (float*)malloc(n * sizeof(float));//______________________________________________________________//void tred2(float* V, float* e, float* d, int n, int ldv);tred2(V_evd, d_evd, e, n, ldv_evd);printf("V_evd-tred2=\n");print_matrix(V_evd, n, n, ldv_evd);printf("\nd-tred2=\n");print_vector(d_evd, n);printf("\ne-tred2=\n");print_vector(e, n);//void tql2(float* V, float* d, float* e, int n, int ldv);tql2(V_evd, d_evd, e, n, ldv_evd);printf("\nV_evd-tql2=\n");print_matrix(V_evd, n, n, ldv_evd);printf("\nd-tql2=\n");print_vector(d_evd, n);printf("\ne-tql2=\n");print_vector(e, n);//_________________________________________________________________________________/**float* d_cevd = nullptr;float2* V_cevd = nullptr;d_cevd = (float*)malloc(N * sizeof(float));V_cevd = (float2*)malloc(ldv_cevd * N * sizeof(float2));*/int ldv_cevd = N;for (int i = 0; i < N; i++)d_cevd[i] = d_evd[2 * i];printf("\nd_cevd=\n");print_vector(d_cevd, N);for (int j = 0; j < N; ++j){int j2 = 2 * j;for (int i = 0; i < N; ++i) {V_cevd[i + j * ldv_cevd].x = V_evd[i + j2 * ldv_evd];//complex<Type>(RV[i][j2], RV[i + N][j2]);V_cevd[i + j * ldv_cevd].y = V_evd[(i + N) + j2 * ldv_evd];}}printf("\nV_cevd=\n");print_matrix_complex(V_cevd, N, N, ldv_cevd);}void print_vector(float* A, int N) {for (int idx = 0; idx < N; idx++) {printf("%7.4f ", A[idx]);}
}
void print_matrix(float* A, int M, int N, 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_matrix_complex(float2* A, int M, int N, int lda) {for (int i = 0; i < M; i++) {for (int j = 0; j < N; j++) {//*(sizeof(complex<float>)/sizeof(float))printf("(%7.4f, %7.4f)", A[i + j * lda].x, A[i + j * lda].y);}printf("\n");}
}

运行：

 

结果分析：

 

matlab:


http://chatgpt.dhexx.cn/article/kessGyPR.shtml

相关文章

Truncated Power Method for Sparse Eigenvalue Problems

Truncated Power Method for Sparse Eigenvalue Problems

目录 算法k的选择\(x\)的初始化代码 抱歉&#xff0c;真的没怎么看懂&#xff0c;当然&#xff0c;估计和我现在没法静下心来好好看也有关系。 算法 想法非常非常简单吧&#xff0c;就是在原来幂法的基础上&#xff0c;每次迭代的时候再加个截断。当然&#xff0c;论文里给出了…
阅读更多...
R语言将向量数据按照行方式转化为矩阵数据（设置参数byrow为TRUE）、计算矩阵数据的特征值（eigenvalue）

R语言将向量数据按照行方式转化为矩阵数据（设置参数byrow为TRUE）、计算矩阵数据的特征值（eigenvalue）

R语言将向量数据按照行方式转化为矩阵数据&#xff08;设置参数byrow为TRUE&#xff09;、计算矩阵数据的特征值&#xff08;eigenvalue&#xff09; 目录 R语言将向量数据按照行方式转化为矩阵数据&#xff08;设置参数byrow为TRUE&#xff09;、计算矩阵数据的特征值&#x…
阅读更多...
Nonlinear Component Analysis as a Kernel Eigenvalue Problem

Nonlinear Component Analysis as a Kernel Eigenvalue Problem

目录 引kernel PCA kernel 的选择性质一些问题代码 Scholkopf B, Smola A J, Muller K, et al. Nonlinear component analysis as a kernel eigenvalue problem[J]. Neural Computation, 1998, 10(5): 1299-1319. 引 普通的PCA将下式进行特征分解&#xff08;用论文的话讲就是对…
阅读更多...
拉格朗日乘数法和KKT条件的直观解释

拉格朗日乘数法和KKT条件的直观解释

拉格朗日乘数法和KKT条件的直观解释 标签&#xff08;空格分隔&#xff09;&#xff1a; 机器学习 linbin 2018-05-10 Abstract 在SVM的推导中&#xff0c;最优化问题是其中的核心&#xff0c;这里我们简单介绍下最优化问题&#xff0c;特别是带有约束的最优化问题&#xff…
阅读更多...
[Math Algorithm] 拉格朗日乘数法

[Math Algorithm] 拉格朗日乘数法

https://www.cnblogs.com/maybe2030/p/4946256.html 阅读目录 1. 拉格朗日乘数法的基本思想2. 数学实例3. 拉格朗日乘数法的基本形态4. 拉格朗日乘数法与KKT条件 拉格朗日乘数法&#xff08;Lagrange Multiplier Method&#xff09;之前听数学老师授课的时候就是一知半解&…
阅读更多...
每天五分钟机器学习算法：拉格朗日乘数法和KKT条件

每天五分钟机器学习算法：拉格朗日乘数法和KKT条件

KKT条件 当我们要求一个函数的极值，同时还有两种类型的约束条件，一种约束条件是等式约束，另外一种约束是不等式约束： x是一个变量（n维，n个样本），我们想要找到使得f(x)最大的x，还要满足上面的约束。此时KKT条件就出来说话了，如果要想让x满足这个条件下的f(x)的最大…
阅读更多...
拉格朗日乘数法及python实现拉格朗日乘数法

拉格朗日乘数法及python实现拉格朗日乘数法

拉格朗日乘数法&#xff08;Lagrange Multiplier Method&#xff09;基本思想 作为一种优化算法&#xff0c;拉格朗日乘子法主要用于解决约束优化问题&#xff0c;它的基本思想就是通过引入拉格朗日乘子来将含有n个变量和k个约束条件的约束优化问题转化为含有&#xff08;nk&am…
阅读更多...
SVM（一）：拉格朗日乘数法详解

SVM（一）：拉格朗日乘数法详解

目录 what直观理解法高数书上的解法 学习SVM的过程中遇到了这个拉格朗日乘数法&#xff0c;之前学高数的时候也学过&#xff0c;不过看到视频里的直观理解法和高数书上的解法有些不同&#xff0c;于是在这里把这两种方法记录下来&#xff0c;也当做是一次理解的过程。 what 先…
阅读更多...
拉格朗日乘数法什么时候考虑端点？解得的点是什么？

拉格朗日乘数法什么时候考虑端点？解得的点是什么？

问题提出 2013年的真题有一道题是用拉格朗日乘数法只能求出来一个点&#xff0c;当时很费解&#xff0c;因此查阅相关资料后&#xff0c;对这部分的知识做一个小总结。 无条件极值和条件极值 首先&#xff0c;在求无条件极值的时候&#xff0c;我们求的是曲面上的极值点。 例…
阅读更多...
拉格朗日乘数法的原理，我用10幅图把它讲清楚

拉格朗日乘数法的原理，我用10幅图把它讲清楚

机器学习是一个目标函数优化问题&#xff0c;给定目标函数f&#xff0c;约束条件会有一般包括以下三类&#xff1a; 仅含等式约束仅含不等式约束等式和不等式约束混合型 当然还有一类没有任何约束条件的最优化问题 关于最优化问题&#xff0c;大都令人比较头疼&#xff0c;首先…
阅读更多...
java 获取随机数方法,java生成随机数的三种方法

java 获取随机数方法,java生成随机数的三种方法

随机数有三种生成方式&#xff1a; 1、通过Math.random()方法 2、通过System.currentTimeMillis()方法获取毫秒数 3、通过Random类 第一种&#xff1a;常用方法Math.random()方法&#xff0c;是获取0-1之间的double类型的小数&#xff0c;在通过int类型墙砖即可 示例&#xff1…
阅读更多...
Java生成随机数的方式

Java生成随机数的方式

目录 Random基础使用优缺点分析 SecureRandom基础使用 总结&#xff1a;持续更新 Random Random 类诞生于 JDK 1.0&#xff0c;它产生的随机数是伪随机数&#xff0c;也就是有规则的随机数。Random 使用的随机算法为 linear congruential pseudorandom number generator (LGC)…
阅读更多...
Java 生成随机数的 5 种方式，你知道几种？

Java 生成随机数的 5 种方式，你知道几种？

点击上方“码农突围”&#xff0c;马上关注 这里是码农充电第一站&#xff0c;回复“666”&#xff0c;获取一份专属大礼包 真爱&#xff0c;请设置“星标”或点个“在看” 作者&#xff1a;专职跑龙套链接&#xff1a;https://www.jianshu.com/p/2f6acd169202 1. Math.random(…
阅读更多...
java 生成随机数 （Random函数）

java 生成随机数 （Random函数）

目录 一、Random是什么&#xff1f; 二、使用步骤 1.引入库 2.创建对象 3.生成随机数 4.完整代码 总结 一、Random是什么&#xff1f; 生成随机数函数 二、使用步骤 1.引入库 代码如下&#xff1a; import java.util.Random; 2.创建对象 代码如下&#xff1a; R…
阅读更多...
Java中生成随机数的4种方式！

Java中生成随机数的4种方式！

在 Java 中&#xff0c;生成随机数的场景有很多&#xff0c;所以本文我们就来盘点一下 4 种生成随机数的方式&#xff0c;以及它们之间的区别和每种生成方式所对应的场景。 1.Random Random 类诞生于 JDK 1.0&#xff0c;它产生的随机数是伪随机数&#xff0c;也就是有规则的…
阅读更多...
谈论SQL注入攻击的重要性

谈论SQL注入攻击的重要性

"SQL注入”是一种利用未过滤/未审核用户输入的攻击方法&#xff08;“缓存溢出”和这个不同&#xff09;&#xff0c;意思就是让应用运行本不应该运行的SQL代码。黑客或者恶搞的用户&#xff0c;利用了程序开发人员在开发的时候没有对SQL进行严格的处理而造成的漏洞&#…
阅读更多...
【SQL注入攻击介绍】

【SQL注入攻击介绍】

目录 前言 本质和危害 分类 注入一般步骤 注入实战 前言 sql注入一直以来都稳居owasp-top10榜首&#xff0c;近年来更是爆出很多的数据库泄露攻击事件&#xff0c;如最近上海某公安存在数据库泄露事件。今天简单的分析以下sql注入的一些特性和方式&#xff1a; owasp-t…
阅读更多...
SQL注入攻击与防护

SQL注入攻击与防护

目录 一、SQL注入攻击概述 1.1 SQL注入概念 1.1.1 标准查询过程 1.1.2 SQL注入定义 1.2 SQL注入根本原因 1.3 SQL注入条件 1.4 SQL注入防范 1.4.1 根本原因&#xff1a;过滤不严 1.4.2 安全设计原则&#xff1a;数据与代码分离 1.5 SQL注入流程 1.6 SQL注入分类 1.…
阅读更多...
使用日志进行调查 - SQL 注入攻击示例

使用日志进行调查 - SQL 注入攻击示例

日志文件是服务器提供的非常有价值的信息。几乎所有服务器、服务和应用程序都提供某种日志记录。日志文件记录在服务或应用程序运行期间发生的事件和操作。 日志文件为我们提供了服务器行为的精确视图以及关键信息&#xff0c;例如何时、如何以及由谁访问服务器。此类信息可以…
阅读更多...
Web—SQL注入攻击

Web—SQL注入攻击

文章目录 一、mysql常用语句二、SQL注入概念1. 产生原因2. 攻击分类 三、攻击流程1. 常用检测语句如何识别SQL注入2. Mysql注入常用函数3. 查询数据的核心语法4. 联合查询5. 报错注入6. 布尔盲注7. 时间盲注8. SQL注入爆库语句9. Sqlmap常用命令 四、常见防护手段及绕过方式1. …
阅读更多...
推荐文章

Copyright @ 2022~2023