Lanzso法を用いて固有ベクトルを求めたいが、うまく行かない
Q&A
Closed
解決したいこと
Lanzso法を用いて固有値、固有ベクトルを求めたいが、うまく行かない。
現状Lanczos法で固有値の計算は最小、第二最小程度まではよく求められているが、そこから固有ベクトルを求めることがうまく行っていない。
以下は参考にしたサイトです。
発生している問題・エラー
Lanczos法でいい固有ベクトルを求める際、プログラムにもあるよう何もLanczosを回して正確な固有ベクトルを求めようとしているが、本来初期ベクトルで固有ベクトルに近いものを選べば収束が早くなるのだが、以下のZを大きくするとどんどん収束から遠くなる（むしろ、Z=０のはじめの周回のほうが固有ベクトルに近いものが出る）。
該当するソースコード
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#define size 800
int M, limit;
void file_write_energy(double Energy[100]){
FILE* fp_energy;
if ((fp_energy = fopen("test_eigenvalue.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
for(int i = 0; i < M; i++){
fprintf(fp_energy, "%lf\n", Energy[i]);
// printf("%lf\n",Energy[i]);
}
fclose(fp_energy);
}
void file_write_vector(double eigenvector_re[size]){
FILE* fp_eigenvector;
if ((fp_eigenvector = fopen("test_eigenvector.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
for(int j = 0; j < size; j++){
fprintf(fp_eigenvector, "%lf\n", eigenvector_re[j]); //基底状態の固有ベクトル
// printf("%lf\n",Energy[i]);
}
fclose(fp_eigenvector);
}
void file_open(double hamiltonian[size][size]){
FILE* fp_hamiltonian;
fp_hamiltonian = fopen("hamiltonian.txt", "r");
if(fp_hamiltonian == NULL) {
printf("ファイルを開くことが出来ませんでした．¥n");
}
for(int l=0; l<size; l++){
for(int i=0;i<size;i++){
fscanf(fp_hamiltonian, "%lf", &(hamiltonian[l][i]) );
}
}
fclose(fp_hamiltonian);
}
//random vector creation
void r_v_c(double vector_q_0[size]){
// 乱数の初期化
srand(time(NULL));
// ランダムな成分からなるベクトルを生成
for (int i = 0; i < size; i++) {
vector_q_0[i] = ((double)rand() / RAND_MAX) * 2.0 - 1.0;
}
double norm = 0.0;
for (int i = 0; i < size; i++) {
norm += vector_q_0[i] * vector_q_0[i];
}
norm = sqrt(norm);
for (int i = 0; i < size; i++) {
vector_q_0[i] /= norm;
}
FILE* fp_rv;
if ((fp_rv = fopen("rv.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
for(int i = 0; i < size; i++){
fprintf(fp_rv, "%lf\n", vector_q_0[i]);
// printf("%lf\n",Energy[i]);
}
fclose(fp_rv);
}
void r_v_c_re(double vector_q_0[size]){
FILE* fp_2;
fp_2 = fopen("rv.txt", "r");
// fp_2 = fopen("test_eigenvector.txt", "r");
if(fp_2 == NULL) {
printf("ファイルを開くことが出来ませんでした．¥n");
}
for(int l = 0; l < size; l++){
fscanf(fp_2, "%lf", &(vector_q_0[l]) );
}
fclose(fp_2);
}
void HV(double hamiltonian[size][size], double vector[size], double out_put_vector[size]){
for(int i = 0; i < size; i++){
for(int j = 0; j < size; j++){
out_put_vector[i] += hamiltonian[i][j] * vector[j];
}
}
}
void matrix(double hamiltonian[size][size], double out_put_vector[size], double vector_q[size]){
HV(hamiltonian, vector_q, out_put_vector);
}
//対格化
void dsyev_(char* JOBZ, char* UPLO, int* n, double* A, int* LDA, double* W, double* WORK, int* LWORK, int* INFO);
void diago(double new_hamiltonian[M][M], double Energy[100], double lanczos_eigenvector[100]){
char JOBZ = 'V';
char UPLO = 'U';
int n = M;
double A[M * M];
for(int i = 0; i < M; i++){
for(int j = 0; j < M; j++){
A[i * M + j] = new_hamiltonian[i][j];
}
}
for(int j = 0; j < M; j++){
lanczos_eigenvector[j] = 0.0;
Energy[j] = 0.0;
}
int LDA = M;
double W[M];
double WORK[6 * M];
int LWORK = 6 * M;
int INFO;
dsyev_(&JOBZ, &UPLO, &n, A, &LDA, W, WORK, &LWORK, &INFO);
for(int i = 0; i < M; i++){
// printf("%lf\n",W[i]);
Energy[i] = W[i];
}
for(int j = 0; j < M; j++){
lanczos_eigenvector[j] = A[0 * M + j]; //基底状態の固有ベクトル
// printf("%d %lf\n", j, lanczos_eigenvector[j]);
}
}
void make_new_hamiltonian(double new_T[M][M], double alpha[100], double beta[100]){
for(int i = 0; i < M; i++){
for(int j = 0; j < M; j++){
new_T[i][j] = 0.0;
}
}
new_T[0][0] = alpha[0], new_T[0][1] = beta[0], new_T[1][0] = beta[0];
for(int i = 0; i < M - 1; i++){
new_T[i][i] = alpha[i], new_T[i][i + 1] = beta[i], new_T[i + 1][i] = beta[i];
}
new_T[M - 1][M - 1] = alpha[M - 1];
}
void kakunin_diago(double alpha[100], double beta[100], double Energy[100],
double lanczos_eigenvector[100], double pre_Energy[2]){
double new_hamiltonian[M][M];
make_new_hamiltonian(new_hamiltonian, alpha, beta);
diago(new_hamiltonian, Energy, lanczos_eigenvector);
}
double break_point_kitei(double Energy[100], double pre_Energy[2]){
double differ = Energy[0] - pre_Energy[0];
pre_Energy[0] = Energy[0];
return differ;
}
double break_point_daiitireiki(double Energy[100], double pre_Energy[2]){
double differ = Energy[1] - pre_Energy[1];
pre_Energy[1] = Energy[1];
return differ;
}
//lanczos eigenvalue
void lanczos_first(double hamiltonian[size][size], double vector_q_0[size], double vector_q_1[size], double a[100], double b[100]){ //k=0
double alpha = 0.0, beta = 0.0;
double Aq[size] = {0.0}, Bq[size] = {0.0};
double bb = 0.0;
matrix(hamiltonian, Aq, vector_q_0);
for(int i = 0; i < size; i++){
alpha += vector_q_0[i] * Aq[i];
}
for(int i = 0; i < size; i++){
Bq[i] = Aq[i] - alpha * vector_q_0[i];
}
for(int i = 0; i < size; i++){
bb += pow(Bq[i], 2);
}
beta = pow(bb, 0.5);
for(int i = 0; i < size; i++){
vector_q_1[i] = Bq[i] / beta;
}
a[0] = alpha, b[0] = beta;
// printf("%lf %lf\n ", alpha, beta);
}
void lanczos_k(double hamiltonian[size][size], double vector_q_0[size], double vector_q_1[size], double vector_q_2[size], double a[100],
double b[100], int lanczos_index){ //0<k<N-2
double alpha = 0.0, beta = 0.0;
double Aq[size] = {0.0}, Bq[size] = {0.0};
double bb = 0.0;
matrix(hamiltonian, Aq, vector_q_1);
for(int i = 0; i < size; i++){
alpha += vector_q_1[i] * Aq[i];
}
for(int i = 0; i < size; i++){
Bq[i] = Aq[i] - b[lanczos_index - 1] * vector_q_0[i] - alpha * vector_q_1[i];
}
for(int i = 0; i < size; i++){
bb += pow(Bq[i], 2);
}
beta = pow(bb, 0.5);
for(int i = 0; i < size; i++){
vector_q_2[i] = Bq[i] / beta;
}
a[lanczos_index] = alpha, b[lanczos_index] = beta;
for(int i = 0; i < size; i++){
vector_q_0[i] = vector_q_1[i];
vector_q_1[i] = vector_q_2[i];
}
// printf("%lf %lf\n ", alpha, beta);
}
void lanczos_step(double hamiltonian[size][size], double vector_q_0[size], double Energy[100], double lanczos_eigenvector[100]){
double static vector_q_1[size] = {0.0}, vector_q_2[size] = {0.0};
double static alpha[100] = {0.0}, beta[100] = {0.0};
double pre_Energy[2] = {0.0}, differ_kitei_value = 10.0, differ_daiitireiki_value = 10.0;
M = 5, limit = 0;
lanczos_first(hamiltonian, vector_q_0, vector_q_1, alpha, beta); //k=0
int lanczos_index = 1, should_break = 0;
while (lanczos_index < 100
&& lanczos_index < size
&& should_break == 0
){
M = lanczos_index;
lanczos_k(hamiltonian, vector_q_0, vector_q_1, vector_q_2, alpha, beta, lanczos_index); //0<k<N-2
(lanczos_index == 5) ? kakunin_diago(alpha, beta, Energy, lanczos_eigenvector, pre_Energy) : 0;
(lanczos_index > 5 && lanczos_index % 2 == 0) ? kakunin_diago(alpha, beta, Energy, lanczos_eigenvector, pre_Energy) : 0;
(lanczos_index > 5 && lanczos_index % 2 == 0) ? (differ_kitei_value = break_point_kitei(Energy, pre_Energy),
differ_daiitireiki_value = break_point_daiitireiki(Energy, pre_Energy)): 0;
should_break = (fabs(differ_kitei_value) < 1e-7 && fabs(differ_daiitireiki_value) < 1e-7);
// printf("%lf %d\n ", lanczos_eigenvector[0], lanczos_index);
lanczos_index++;
}
limit = M;
}
void make_eigenvector(double eigenvector_re[size], double vector_q_1[size], double lanczos_eigenvector[100], int index){
for(int i = 0; i < size; i++){
eigenvector_re[i] += lanczos_eigenvector[index] * vector_q_1[i];
}
}
//lanczos eigen vector
void lanczos_first_re(double hamiltonian[size][size], double vector_q_0[size], double vector_q_1[size], double a[100], double b[100],
double eigenvector_re[size], double lanczos_eigenvector[100]){ //k=0
double alpha = 0.0, beta = 0.0;
double Aq[size] = {0.0}, Bq[size] = {0.0};
double bb = 0.0;
matrix(hamiltonian, Aq, vector_q_0);
for(int i = 0; i < size; i++){
alpha += vector_q_0[i] * Aq[i];
}
for(int i = 0; i < size; i++){
Bq[i] = Aq[i] - alpha * vector_q_0[i];
}
for(int i = 0; i < size; i++){
bb += pow(Bq[i], 2);
}
beta = pow(bb, 0.5);
for(int i = 0; i < size; i++){
vector_q_1[i] = Bq[i] / beta;
}
a[0] = alpha, b[0] = beta;
make_eigenvector(eigenvector_re, vector_q_1, lanczos_eigenvector, 0);
}
void lanczos_k_re(double hamiltonian[size][size], double vector_q_0[size], double vector_q_1[size], double vector_q_2[size], double a[100], double b[100],
double eigenvector_re[size], double lanczos_eigenvector[100], int index){ //0<k<N-2
double alpha = 0.0, beta = 0.0;
double Aq[size] = {0.0}, Bq[size] = {0.0};
double bb = 0.0;
matrix(hamiltonian, Aq, vector_q_1);
for(int i = 0; i < size; i++){
alpha += vector_q_1[i] * Aq[i];
}
for(int i = 0; i < size; i++){
Bq[i] = Aq[i] - b[index - 1] * vector_q_0[i] - alpha * vector_q_1[i];
}
for(int i = 0; i < size; i++){
bb += pow(Bq[i], 2);
}
beta = pow(bb, 0.5);
for(int i = 0; i < size; i++){
vector_q_2[i] = Bq[i] / beta;
}
a[index] = alpha, b[index] = beta;
for(int i = 0; i < size; i++){
vector_q_0[i] = vector_q_1[i];
vector_q_1[i] = vector_q_2[i];
}
make_eigenvector(eigenvector_re, vector_q_1, lanczos_eigenvector, index);
}
void lanczos_step_re(double hamiltonian[size][size], double vector_q_0[size], double eigenvector_re[size], double lanczos_eigenvector[100]){
double vector_q_1[size] = {0.0}, vector_q_2[size] = {0.0};
double alpha[100] = {0.0}, beta[100] = {0.0};
lanczos_first_re(hamiltonian, vector_q_0, vector_q_1, alpha, beta, eigenvector_re, lanczos_eigenvector); //k=0
int index = 1;
while (index <= limit){
lanczos_k_re(hamiltonian, vector_q_0, vector_q_1, vector_q_2, alpha, beta, eigenvector_re, lanczos_eigenvector, index); //0<k<N-2
index++;
}
// printf("%d\n", limit);
}
void normalization(double vector[size]){
double norm = 0.0;
for(int i = 0; i < size; i++){
norm += pow(vector[i], 2);
}
norm = pow(norm, 0.5);
for(int i = 0; i < size; i++){
vector[i] = vector[i] / norm;
}
}
int main(){
double vector_q_0[size] = {0.0}, Energy[100] = {0.0}, eigenvector_re[size] = {0.0};
double hamiltonian[size][size] = {0.0};
file_open(hamiltonian);
// r_v_c(vector_q_0); //random vector_q_0
r_v_c_re(vector_q_0); //saiki vector_q_0
int z = 0;
while(z < 5){ //←この更新数が多くなると不確かになる
double vector_Q[size] = {0.0}, lanczos_eigenvector[100] = {0.0};
M = 5, limit = 0;
for(int i = 0; i < size; i++){
vector_Q[i] = vector_q_0[i];
}
lanczos_step(hamiltonian, vector_q_0, Energy, lanczos_eigenvector);
for(int i = 0; i < size; i++){
vector_q_0[i] = vector_Q[i];
}
lanczos_step_re(hamiltonian, vector_q_0, eigenvector_re, lanczos_eigenvector);
for(int i = 0; i < size; i++){
vector_q_0[i] = eigenvector_re[i];
eigenvector_re[i] = 0.0;
}
normalization(vector_q_0);
z++;
}
file_write_energy(Energy);
file_write_vector(vector_q_0);
return 0;
}
自分で試したこと
まず、二回のLanczos法で同じベクトルや固有値が計算できているかを調べたが、全く同じであることが確認できた。数値の受け渡しも特に問題がないことがわかっている。
固有値に関してはLanczos法の反復が途中で終わることから１００回までの間で
きちんと1e-7の精度で収束しているので上のlanczos_stepは問題ないと考えている。では下じゃないか、と思ったが、いろいろ変えてみても何も良くならなかったので、他人のお力をお借りしたいです。
環境はよくわかりませんが、ubuntuでうｇ
0
