500分急急求算矩阵特征值的改良HQR算法

[goby]

现在求矩阵特征值一般都用HQR算法，但好像通常能找到的HQR算法都有个问题：就是若矩阵类似正交矩阵时算法失效。不知道有没有改良的算法，最好能给源代码。
mathematica和matlab的HQR算法好像就没问题。
（由于最多只能一次给100分，需要开几个帖子才能给，请见谅）

[saint001]

没见过HQR算法
我以前写的一个求矩阵特征值的程序是先化成Hessenburg阵，然后用QR
不知道什么是HQR

[saint001]

这是我从Numerical Recipes in C上找到的hqr函数，稍作修改，运行应该是正确的
#include "math.h"
#include "stdio.h"

#define SIGN(a,b) ((b) > 0 ? fabs(a) : -fabs(a))

void hqr(double **a,int n,double *wr,double *wi)
{
int nn,m,l,k,j,its,i,mmin;
double z,y,x,w,v,u,t,s,r,q,p,anorm;

anorm=fabs(a[1][1]);
for (i=2;i<=n;i++)
for (j=(i-1);j<=n;j++)
anorm += fabs(a[i][j]);
nn=n;
t=0.0;
while (nn >= 1) {
its=0;
do {
for (l=nn;l>=2;l--) {
s=fabs(a[l-1][l-1])+fabs(a[l][l]);
if (s == 0.0) s=anorm;
if (fabs(a[l][l-1]) + s == s) break;
}
x=a[nn][nn];
if (l == nn) {
wr[nn]=x+t;
wi[nn--]=0.0;
} else {
y=a[nn-1][nn-1];
w=a[nn][nn-1]*a[nn-1][nn];
if (l == (nn-1)) {
p=0.5*(y-x);
q=p*p+w;
z=sqrt(fabs(q));
x += t;
if (q >= 0.0) {
z=p+SIGN(z,p);
wr[nn-1]=wr[nn]=x+z;
if (z) wr[nn]=x-w/z;
wi[nn-1]=wi[nn]=0.0;
} else {
wr[nn-1]=wr[nn]=x+p;
wi[nn-1]= -(wi[nn]=z);
}
nn -= 2;
} else {
if (its == 30) printf("Too many iterations in HQR");
if (its == 10 || its == 20) {
t += x;
for (i=1;i<=nn;i++) a[i][i] -= x;
s=fabs(a[nn][nn-1])+fabs(a[nn-1][nn-2]);
y=x=0.75*s;
w = -0.4375*s*s;
}
++its;
for (m=(nn-2);m>=l;m--) {
z=a[m][m];
r=x-z;
s=y-z;
p=(r*s-w)/a[m+1][m]+a[m][m+1];
q=a[m+1][m+1]-z-r-s;
r=a[m+2][m+1];
s=fabs(p)+fabs(q)+fabs(r);
p /= s;
q /= s;
r /= s;
if (m == l) break;
u=fabs(a[m][m-1])*(fabs(q)+fabs(r));
v=fabs(p)*(fabs(a[m-1][m-1])+fabs(z)+fabs(a[m+1][m+1]));
if (u+v == v) break;
}
for (i=m+2;i<=nn;i++) {
a[i][i-2]=0.0;
if (i != (m+2)) a[i][i-3]=0.0;
}
for (k=m;k<=nn-1;k++) {
if (k != m) {
p=a[k][k-1];
q=a[k+1][k-1];
r=0.0;
if (k != (nn-1)) r=a[k+2][k-1];
if (x=fabs(p)+fabs(q)+fabs(r)) {
p /= x;
q /= x;
r /= x;
}
}
if (s=SIGN(sqrt(p*p+q*q+r*r),p)) {
if (k == m) {
if (l != m)
a[k][k-1] = -a[k][k-1];
} else
a[k][k-1] = -s*x;
p += s;
x=p/s;
y=q/s;
z=r/s;
q /= p;
r /= p;
for (j=k;j<=nn;j++) {
p=a[k][j]+q*a[k+1][j];
if (k != (nn-1)) {
p += r*a[k+2][j];
a[k+2][j] -= p*z;
}
a[k+1][j] -= p*y;
a[k][j] -= p*x;
}
mmin = nn<k+3 ? nn : k+3;
for (i=l;i<=mmin;i++) {
p=x*a[i][k]+y*a[i][k+1];
if (k != (nn-1)) {
p += z*a[i][k+2];
a[i][k+2] -= p*r;
}
a[i][k+1] -= p*q;
a[i][k] -= p;
}
}
}
}
}
} while (l < nn-1);
}
}


int main(int argc, char* argv[])
{
int i,j,n;
double **a,*wr,*wi;
FILE *fp=fopen("d:\\hqr.txt","r");
fscanf(fp,"%d",&n);
a=new double*[n+1];
wr=new double[n+1];
wi=new double[n+1];
for(i=0;i<n+1;i++)
a[i]=new double[n+1];
for(i=1;i<n+1;i++)
for(j=1;j<n+1;j++)
fscanf(fp,"%lf",a[i]+j);
fclose(fp);

printf("原矩阵\n");
for(i=1;i<n+1;i++)
{
for(j=1;j<n+1;j++)
printf("%f\t",a[i][j]);
printf("\n");
}
printf("\n");
hqr(a,n,wr,wi);

printf("Hessenburg阵\n");
for(i=1;i<n+1;i++)
{
for(j=1;j<n+1;j++)
printf("%f\t",a[i][j]);
printf("\n");
}
printf("\n");

printf("特征值\n");
for(i=1;i<n+1;i++)
printf("%f\t+i*%f\n",wr[i],wi[i]);
return 0;
}

d:/hqr.txt中的数据
3
-12 3 3
3 1 2
3 -2 7

[saint001]

就是有问题
不是正交矩阵的时候运算的精度非常高，和matlab求的特征值一样
但取A= -0.40824829046386 0.70710678118655 -0.57735026918963
0.81649658092773 0 -0.57735026918963
-0.40824829046386 -0.70710678118655 -0.57735026918963
时(一个旋转矩阵)程序求出的结果是
Hessenburg阵
0.894439 -0.257957 0.038804
0.000000 -0.831438 -0.388133
0.000000 0.154166 -1.048599

特征值
0.894439 +i*0.000000
-0.940019 +i*-0.219197
-0.940019 +i*0.219197
matlab求的结果是
1.00000000000000
-0.99279927982674 + 0.11978977408568i
-0.99279927982674 - 0.11978977408568i
但如果把a的每一项都加上1
0.59175170953614 1.70710678118655 0.42264973081037
1.81649658092773 1.00000000000000 0.42264973081037
0.59175170953614 0.29289321881345 0.42264973081037
结果就一样了

[merelynobody]

解决了
换个名字up一下

[saint001]

用HQR方法求解特征值
第一步均衡，第二步化为Ｈｅｓｓｅｎｂｕｒｇ阵，第三步HQR
从Numberical Recipes in C有说明，三个函数，我去年编的就是这个

#include "math.h"
#include "stdio.h"

#define SIGN(a,b) ((b) > 0 ? fabs(a) : -fabs(a))
#define SWAP(g,h) {y=(g);(g)=(h);(h)=y;}
#define RADIX 2.0

void balanc(double **a, int n)
{
int last,j,i;
double s,r,g,f,c,sqrdx;
sqrdx=RADIX*RADIX;
last=0;
while (last == 0)
{
last=1;
for (i=1;i<=n;i++)
{
r=c=0.0;
for (j=1;j<=n;j++)
if (j != i)
{
c += fabs(a[j][i]);
r += fabs(a[i][j]);
}
if (c && r)
{
g=r/RADIX;
f=1.0;
s=c+r;
while (c<g)
{
f *= RADIX;
c *= sqrdx;
}
g=r*RADIX;
while (c>g)
{
f /= RADIX;
c /= sqrdx;
}
if ((c+r)/f < 0.95*s)
{
last=0;
g=1.0/f;
for (j=1;j<=n;j++)
a[i][j] *= g;
for (j=1;j<=n;j++)
a[j][i] *= f;
}
}
}
}
}

void elmhes(double **a,int n)
{
int m,j,i;
double y,x;

for (m=2;m<n;m++) {
x=0.0;
i=m;
for (j=m;j<=n;j++) {
if (fabs(a[j][m-1]) > fabs(x)) {
x=a[j][m-1];
i=j;
}
}
if (i != m) {
for (j=m-1;j<=n;j++) SWAP(a[i][j],a[m][j])
for (j=1;j<=n;j++) SWAP(a[j][i],a[j][m])
}
if (x) {
for (i=m+1;i<=n;i++) {
if (y=a[i][m-1]) {
y /= x;
a[i][m-1]=y;
for (j=m;j<=n;j++)
a[i][j] -= y*a[m][j];
for (j=1;j<=n;j++)
a[j][m] += y*a[j][i];
}
}
}
}
}

void hqr(double **a,int n,double *wr,double *wi)
{
int nn,m,l,k,j,its,i,mmin;
double z,y,x,w,v,u,t,s,r,q,p,anorm;

anorm=fabs(a[1][1]);
for (i=2;i<=n;i++)
for (j=(i-1);j<=n;j++)
anorm += fabs(a[i][j]);
nn=n;
t=0.0;
while (nn >= 1) {
its=0;
do {
for (l=nn;l>=2;l--) {
s=fabs(a[l-1][l-1])+fabs(a[l][l]);
if (s == 0.0) s=anorm;
if (fabs(a[l][l-1]) + s == s) break;
}
x=a[nn][nn];
if (l == nn) {
wr[nn]=x+t;
wi[nn--]=0.0;
} else {
y=a[nn-1][nn-1];
w=a[nn][nn-1]*a[nn-1][nn];
if (l == (nn-1)) {
p=0.5*(y-x);
q=p*p+w;
z=sqrt(fabs(q));
x += t;
if (q >= 0.0) {
z=p+SIGN(z,p);
wr[nn-1]=wr[nn]=x+z;
if (z) wr[nn]=x-w/z;
wi[nn-1]=wi[nn]=0.0;
} else {
wr[nn-1]=wr[nn]=x+p;
wi[nn-1]= -(wi[nn]=z);
}
nn -= 2;
} else {
if (its == 30) printf("Too many iterations in HQR");
if (its == 10 || its == 20) {
t += x;
for (i=1;i<=nn;i++) a[i][i] -= x;
s=fabs(a[nn][nn-1])+fabs(a[nn-1][nn-2]);
y=x=0.75*s;
w = -0.4375*s*s;
}
++its;
for (m=(nn-2);m>=l;m--) {
z=a[m][m];
r=x-z;
s=y-z;
p=(r*s-w)/a[m+1][m]+a[m][m+1];
q=a[m+1][m+1]-z-r-s;
r=a[m+2][m+1];
s=fabs(p)+fabs(q)+fabs(r);
p /= s;
q /= s;
r /= s;
if (m == l) break;
u=fabs(a[m][m-1])*(fabs(q)+fabs(r));
v=fabs(p)*(fabs(a[m-1][m-1])+fabs(z)+fabs(a[m+1][m+1]));
if (u+v == v) break;
}
for (i=m+2;i<=nn;i++) {
a[i][i-2]=0.0;
if (i != (m+2)) a[i][i-3]=0.0;
}
for (k=m;k<=nn-1;k++) {
if (k != m) {
p=a[k][k-1];
q=a[k+1][k-1];
r=0.0;
if (k != (nn-1)) r=a[k+2][k-1];
if (x=fabs(p)+fabs(q)+fabs(r)) {
p /= x;
q /= x;
r /= x;
}
}
if (s=SIGN(sqrt(p*p+q*q+r*r),p)) {
if (k == m) {
if (l != m)
a[k][k-1] = -a[k][k-1];
} else
a[k][k-1] = -s*x;
p += s;
x=p/s;
y=q/s;
z=r/s;
q /= p;
r /= p;
for (j=k;j<=nn;j++) {
p=a[k][j]+q*a[k+1][j];
if (k != (nn-1)) {
p += r*a[k+2][j];
a[k+2][j] -= p*z;
}
a[k+1][j] -= p*y;
a[k][j] -= p*x;
}
mmin = nn<k+3 ? nn : k+3;
for (i=l;i<=mmin;i++) {
p=x*a[i][k]+y*a[i][k+1];
if (k != (nn-1)) {
p += z*a[i][k+2];
a[i][k+2] -= p*r;
}
a[i][k+1] -= p*q;
a[i][k] -= p;
}
}
}
}
}
} while (l < nn-1);
}
}


int main(int argc, char* argv[])
{
int i,j,n;
double **a,*wr,*wi;
FILE *fp=fopen("d:\\hqr.txt","r");
fscanf(fp,"%d",&n);
a=new double*[n+1];
wr=new double[n+1];
wi=new double[n+1];
for(i=0;i<n+1;i++)
a[i]=new double[n+1];
for(i=1;i<n+1;i++)
for(j=1;j<n+1;j++)
fscanf(fp,"%lf",a[i]+j);
fclose(fp);

printf("原矩阵\n");
for(i=1;i<n+1;i++)
{
for(j=1;j<n+1;j++)
printf("%f\t",a[i][j]);
printf("\n");
}
printf("\n");
//关键的三个函数
balanc(a,n);
elmhes(a,n);
hqr(a,n,wr,wi);

printf("特征值\n");
for(i=1;i<n+1;i++)
printf("%f\t+i*%f\n",wr[i],wi[i]);
return 0;
}

[saint001]

现在对正交矩阵
-0.40824829046386 0.70710678118655 -0.57735026918963
0.81649658092773 0 -0.57735026918963
-0.40824829046386 -0.70710678118655 -0.57735026918963
求解的结果就是一样的了
可能关键是balanc函数

[saint001]

可是不用balanc函数也得出正确的结果，虽然可能迭代次数多了一些
Numerical Recipes in C的pdf文档可以从网上下载
http://www.library.cornell.edu/nr/bookcpdf.html
但是网站上不提供程序(其实下载的文档里有大部分的程序，balanc函数就是从书上拷下来的)
我这里有以前从matlab大观园下载的代码，好像现在不能下了，如果想看的话可以发给你

[goby]

好的，谢谢你
如果能提供代码更加不胜感激：[e-mail]
另外，我需要另开几个空贴以便继续给分

[saint001]

已经发过去了