lanczos finished, kryleig update

LectureCodes/Evp/kryleig/CMakeLists.txt

+add_subdirectory(Eigen)

LectureCodes/Evp/kryleig/Eigen/CMakeLists.txt

+project(kryleig)
+cmake_minimum_required(VERSION 3.0.2)
+
+include_directories(../../utils/)
+
+add_executable_numcse(main main.cpp)

LectureCodes/Evp/kryleig/Eigen/kryleig.hpp

+#include <Eigen/Dense>
+#include <Eigen/Eigenvalues>
+
+
+// Ritz projection onto Krylov subspace. An orthonormal basis of \Blue{$\Kryl{m}{\VA}{\mathbf{1}}$} is assembled into the columns of \Blue{$\VV$}. 
+void kryleig(const Eigen::MatrixXd &A, int m, Eigen::MatrixXd &V, Eigen::VectorXd &d)
+{
+	int n = A.rows();
+	V.resize(n,m);
+	d.resize(m);
+   	V.col(0) = Eigen::VectorXd::LinSpaced(n,1,n);
+	V.col(0).normalize();
+
+	for (int l=1; l<m; ++l)
+	{
+		V.col(l) = A*V.col(l-1);
+
+		auto qr = V.householderQr();
+		Eigen::MatrixXd Q = qr.householderQ() * Eigen::MatrixXd::Identity(n,l+1); // economy size qr decomposition
+
+		Eigen::EigenSolver<Eigen::MatrixXd> esolver(Q.transpose()*A*Q);
+
+		V.leftCols(l+1) = Q*esolver.eigenvectors().real();
+		d.head(l+1) = esolver.eigenvalues().real();
+	}
+}

LectureCodes/Evp/kryleig/Eigen/main.cpp

+#include <iostream>
+#include <Eigen/Dense>
+#include <figure.hpp>
+#include "eigensolversort.hpp"
+#include "kryleig.hpp"
+#include "spdiags.hpp"
+
+
+int main()
+{
+	// Generate matrix
+	
+	int n = 100;
+	Eigen::MatrixXd M = Eigen::MatrixXd::Zero(n,n);
+	M.diagonal(-1) = -0.5*Eigen::VectorXd::Ones(n-1);
+	M.diagonal(0) =   2.*Eigen::VectorXd::Ones(n);
+	M.diagonal(1) = -1.5*Eigen::VectorXd::Ones(n-1);
+
+	auto qr = M.householderQr();   // generate orthogonal matrix
+	Eigen::MatrixXd Q = qr.householderQ();
+
+	// eigenvalues \Blue{$1,2,3,\ldots,100$}
+	Eigen::VectorXd diag= Eigen::VectorXd::LinSpaced(n,1,n);
+	Eigen::MatrixXd A = Q.transpose()*diag.asDiagonal()*Q;
+
+	Eigen::EigenSolver<Eigen::MatrixXd> esolver(A);
+
+	Eigen::VectorXd d0; // eigenvalues
+	Eigen::MatrixXd V0; // eigenvectors
+	std::tie(d0, V0) = eigensolversort(esolver);
+
+	Eigen::VectorXd ev0(6);
+	ev0.head(3) = d0.head(3);
+	ev0.tail(3) = d0.tail(3);
+
+	int N = 30-5;
+	Eigen::MatrixXd errev(N,13);
+
+	// kryleig return values 
+	Eigen::VectorXd d;
+	Eigen::MatrixXd V;
+	Eigen::VectorXd ev(6);
+
+	for (int i=0; i<N; ++i)
+	{
+		const int m = i+6;
+		std::cout << "subspace dimension m = " << m << std::endl;
+		kryleig(A,m, V, d);
+		std::sort(d.data(),d.data()+d.size());	   
+		ev.head(3) = d.head(3);
+		ev.tail(3) = d.tail(3);
+
+		errev(i,0) = m;
+		errev.block(i,1,1,6) = ev.transpose();
+		errev.block(i,7,1,6) = (ev-ev0).cwiseAbs().transpose();
+	}
+
+	{
+		mgl::Figure fig;
+		fig.plot(errev.col(0), errev.col(1), "r+");
+		fig.plot(errev.col(0), errev.col(2), "m+");
+		fig.plot(errev.col(0), errev.col(3), "b+");
+		fig.title("Krylov subspace iteration - smallest ritz values");
+		fig.xlabel("Dimension m of Krylov space");
+		fig.ylabel("Ritz value");
+		fig.setFontSize(3);
+		fig.save("kryleigmin");
+	}
+
+	{
+		mgl::Figure fig;
+		fig.plot(errev.col(0), errev.col(6), "r^");
+		fig.plot(errev.col(0), errev.col(5), "m^");
+		fig.plot(errev.col(0), errev.col(4), "b^");
+		fig.title("Krylov subspace iteration - largest ritz values");
+		fig.xlabel("Dimension m of Krylov space");
+		fig.ylabel("Ritz value");
+		fig.setFontSize(3);
+		fig.save("kryleigmax");
+	}
+
+	{
+		mgl::Figure fig;
+		fig.setlog(false, true);
+		fig.plot(errev.col(0), errev.col(7), "r+");
+		fig.plot(errev.col(0), errev.col(8), "m+");
+		fig.plot(errev.col(0), errev.col(9), "b+");
+		fig.title("Krylov subspace iteration - error of smallest ritz values");
+		fig.xlabel("Dimension m of Krylov space");
+		fig.ylabel("Error of Ritz value");
+		fig.setFontSize(3);
+		fig.save("kryleigminerr");
+	}
+
+	{
+		mgl::Figure fig;
+		fig.setlog(false, true);
+		fig.plot(errev.col(0), errev.col(10), "r^");
+		fig.plot(errev.col(0), errev.col(11), "m^");
+		fig.plot(errev.col(0), errev.col(12), "b^");
+		fig.title("Krylov subspace iteration - error of largest ritz values");
+		fig.xlabel("Dimension m of Krylov space");
+		fig.ylabel("Error of Ritz value");
+		fig.setFontSize(3);
+		fig.save("kryleigmaxerr");
+	}
+
+}

LectureCodes/Evp/kryleig/Matlab/kryleig.m

+function [V,D] = kryleig(A,m)
+% Ritz projection onto Krylov subspace. An orthonormal basis of \Blue{$\Kryl{m}{\VA}{\mathbf{1}}$} is assembled into the columns of \Blue{$\VV$}. 
+n = size(A,1); V = (1:n)'; V = V/norm(V); 
+for l=1:m-1
+  V = [V,A*V(:,end)]; [Q,R] = qr(V,0); 
+  [W,D] = eig(Q'*A*Q); V = Q*W;
+end

LectureCodes/Evp/kryleig/Matlab/kryleigd.m

+% MATLAB script for ex:kryleig
+n=100; M=gallery('tridiag',-0.5*ones(n-1,1),2*ones(n,1),-1.5*ones(n-1,1));
+[Q,R]=qr(M); A=Q'*diag(1:n)*Q; % eigenvalues \Blue{$1,2,3,\ldots,100$}
+
+[V0,D0] = eig(A); d0 = sort(diag(D0)); 
+ev0 = [d0(1:3);d0(end-2:end)];
+
+errev = [];
+for m=6:30
+  [V,D] = kryleig(A,m); d = sort(diag(D));
+  ev =  [d(1:3);d(end-2:end)];
+  errev = [errev; m , ev', abs((ev-ev0)')];
+end
+
+figure; 
+plot(errev(:,1),errev(:,2),'r+',errev(:,1),errev(:,3),'m+',errev(:,1),errev(:,4),'b+');
+axis([5 31 -1 41]);
+xlabel('{\bf dimension m of Krylov space}','fontsize',14);
+ylabel('{\bf Ritz value}','fontsize',14);
+legend('\mu_1','\mu_2','\mu_3','location','northeast');
+
+print -depsc2 '../PICTURES/kryleigmin.eps';
+
+figure;
+plot(errev(:,1),errev(:,7),'r^',errev(:,1),errev(:,6),'m^',errev(:,1),errev(:,5),'b^');
+axis([5 31 70 101]);
+xlabel('{\bf dimension m of Krylov space}','fontsize',14);
+ylabel('{\bf Ritz value}','fontsize',14);
+legend('\mu_{m}','\mu_{m-1}','\mu_{m-2}','location','southeast');
+
+print -depsc2 '../PICTURES/kryleigmax.eps';
+
+
+figure;
+semilogy(errev(:,1),errev(:,8),'r+',errev(:,1),errev(:,9),'m+',errev(:,1),errev(:,10),'b+');
+ax = axis; axis([5 31 ax(3) ax(4)]);
+xlabel('{\bf dimension m of Krylov space}','fontsize',14);
+ylabel('{\bf error of Ritz value}','fontsize',14);
+legend('|\lambda_1-\mu_1|','|\lambda_2-\mu_2|','|\lambda_2-\mu_3|',...
+       'location','northeast');
+
+print -depsc2 '../PICTURES/kryleigminerr.eps';
+
+figure;
+semilogy(errev(:,1),errev(:,13),'r^',errev(:,1),errev(:,12),'m^',errev(:,1),errev(:,11),'b^');
+ax = axis; axis([5 31 ax(3) ax(4)]);
+xlabel('{\bf dimension m of Krylov space}','fontsize',14);
+ylabel('{\bf Ritz value}','fontsize',14);
+legend('|\lambda_m-\mu_{m}|','|\lambda_{m-1}-\mu_{m-1}|','|\lambda_{m-2}-\mu_{m-2}|','location','northeast');
+
+print -depsc2 '../PICTURES/kryleigmaxerr.eps';
+
+
+
+

LectureCodes/Evp/lanczos/CMakeLists.txt

+add_subdirectory(Eigen)

LectureCodes/Evp/lanczos/Eigen/CMakeLists.txt

+project(lanczos)
+cmake_minimum_required(VERSION 2.8)
+
+include_directories(../../utils/)
+add_executable_numcse(main main.cpp)

LectureCodes/Evp/lanczos/Eigen/lanczos.hpp

+#include <Eigen/Dense>
+
+/* Note: this implementation of the Lanczos process also records the orthonormal CG residuals in the columns of the matrix \texttt{V}, which is not needed when only eigenvalue approximations are desired. */
+void lanczos(const Eigen::MatrixXd& A, int k, const Eigen::VectorXd &z0, Eigen::MatrixXd &V, Eigen::VectorXd& alpha, Eigen::VectorXd& beta)
+{
+	V.resize(A.rows(),k+1);
+	alpha.resize(k);
+	beta.resize(k);
+
+	V.col(0) = z0 / z0.norm();
+
+	// Vectors storing entries of tridiagonal matrix \eqref{eq:tdcg}
+	alpha = Eigen::VectorXd::Zero(k,1); 
+	beta = Eigen::VectorXd::Zero(k,1);
+
+	for (int j=0; j<k; ++j)
+	{
+		Eigen::VectorXd q = A*V.col(j);
+		alpha(j) = q.dot(V.col(j));
+		Eigen::VectorXd w = q - alpha(j)*V.col(j);
+
+		if (j > 0) w = w - beta(j-1)*V.col(j-1);
+
+		beta(j) = w.norm();
+		V.col(j+1) = w/beta(j);
+	}
+	beta = beta.head(k-1);
+}

LectureCodes/Evp/lanczos/Eigen/lanczosev.hpp

+#include <Eigen/Dense>
+#include <Eigen/Eigenvalues>
+#include <eigensolversort.hpp>
+
+Eigen::MatrixXd lanczosev(const Eigen::MatrixXd& A, int k, int j, const Eigen::VectorXd& w)
+{
+	Eigen::MatrixXd V(A.rows(), k+1);
+	V.col(0) = w / w.norm();
+
+	Eigen::VectorXd alpha = Eigen::VectorXd::Zero(k);
+	Eigen::VectorXd beta = Eigen::VectorXd::Zero(k);
+
+	Eigen::MatrixXd result(k, j+1);
+	Eigen::VectorXd vt;
+
+	for (int l=0; l<k; ++l)
+	{
+		vt = A*V.col(l);
+		if (l > 0) vt = vt - beta(l-1)*V.col(l-1);
+
+		alpha(l) = vt.dot(V.col(l));
+		vt = vt - alpha(l)*V.col(l);
+		beta(l) = vt.norm();
+		V.col(l+1) = vt/beta(l);
+
+		Eigen::MatrixXd T2 = Eigen::MatrixXd::Zero(l+1,l+1);
+
+		// create diagonals
+		if (l>1)
+		{
+			Eigen::VectorXd dl(l);
+			Eigen::VectorXd du(l);
+
+			du = beta.head(l);
+			dl = beta.head(l);
+
+			T2.diagonal(-1) = dl;
+			T2.diagonal(1) = du;
+		}
+
+		T2.diagonal(0) = alpha.head(l+1);
+
+		// eigendecomposition
+		Eigen::EigenSolver<Eigen::MatrixXd> esolver(T2);
+		Eigen::MatrixXd U; 
+		Eigen::VectorXd ev; 
+		std::tie(ev, U) = eigensolversort(esolver);
+
+		Eigen::VectorXd ev2(j);
+
+		if (j > l)
+		{
+			ev2.head(j-l-1) = Eigen::VectorXd::Zero(j-l-1);
+			ev2.tail(l+1) = ev;
+		}
+		else
+		{
+			ev2 = ev.tail(j);
+		}
+		result(l,0) = l;
+		result.row(l).tail(j) = ev2;
+	}
+
+	return result;
+}

LectureCodes/Evp/lanczos/Eigen/main.cpp

+#include <iostream>
+#include <Eigen/Dense>
+#include <figure.hpp>
+#include "lanczos.hpp"
+#include "lanczosev.hpp"
+
+#include <Eigen/Eigenvalues>
+#include <eigensolversort.hpp>
+
+int main()
+{
+	// create test matrix
+	
+	int n = 100;
+	Eigen::MatrixXd M = Eigen::MatrixXd::Zero(n,n);
+	M.diagonal(-1) = -0.5*Eigen::VectorXd::Ones(n-1);
+	M.diagonal(0) =   2.*Eigen::VectorXd::Ones(n);
+	M.diagonal(1) = -1.5*Eigen::VectorXd::Ones(n-1);
+
+	auto qr = M.householderQr();   // generate orthogonal matrix
+	Eigen::MatrixXd Q = qr.householderQ();
+
+	// eigenvalues \Blue{$1,2,3,\ldots,100$}
+	Eigen::VectorXd diag= Eigen::VectorXd::LinSpaced(n,1,n);
+	Eigen::MatrixXd A = Q.transpose()*diag.asDiagonal()*Q;
+
+	// test lancos
+	int k = 30;
+	Eigen::VectorXd z0 = Eigen::VectorXd::Ones(n);
+	Eigen::VectorXd alpha;
+	Eigen::VectorXd beta;
+	Eigen::MatrixXd V;
+	lanczos(A, k, z0, V, alpha, beta);
+
+	std::cout << "alpha:" << std::endl;
+	std::cout << alpha << std::endl;
+	std::cout << "beta:" << std::endl;
+	std::cout << beta << std::endl;
+
+
+	// test lancosev 
+	Eigen::EigenSolver<Eigen::MatrixXd> esolver(A);
+	Eigen::MatrixXd U; 
+	Eigen::VectorXd ev; 
+	std::tie(ev, U) = eigensolversort(esolver);
+
+	
+	int j = 5;
+	Eigen::VectorXd ones = Eigen::VectorXd::Ones(100);
+	Eigen::MatrixXd res = lanczosev(A,k,j,ones);
+
+	{
+		mgl::Figure fig;
+		Eigen::VectorXd y1 = res.col(j) - Eigen::MatrixXd::Constant(k, 1, ev(n-1));
+		Eigen::VectorXd y2 = res.col(j-1) - Eigen::MatrixXd::Constant(k, 1, ev(n-2));
+		Eigen::VectorXd y3 = res.col(j-2) - Eigen::MatrixXd::Constant(k, 1, ev(n-3));
+		Eigen::VectorXd y4 = res.col(j-3) - Eigen::MatrixXd::Constant(k, 1, ev(n-4));
+		Eigen::VectorXd x = Eigen::VectorXd::LinSpaced(k,1,k);
+		fig.setlog(false, true);
+		fig.plot(x, 		  y1.cwiseAbs(), 		   "r+").label("lamda n");
+		fig.plot(x.tail(k-1), y2.cwiseAbs().tail(k-1), "m-^").label("lamda n-1");
+		fig.plot(x.tail(k-2), y3.cwiseAbs().tail(k-2), "b-o").label("lamda n-2");
+		fig.plot(x.tail(k-3), y4.cwiseAbs().tail(k-3), "k-*").label("lamda n-3");
+		
+		fig.legend(1,1);
+		fig.xlabel("step of Lanzcos process");
+		fig.ylabel("|Ritz value-eigenvalue|");
+		fig.setFontSize(3);
+		fig.save("lanczosev1");
+	}
+}

LectureCodes/Evp/lanczos/Matlab/lanczos.m

+function [V,alph,bet] = lanczos(A,k,z0)
+% Note: this implementation of the Lanczos process also records the orthonormal CG residuals in the columns of the matrix \texttt{V}, which is not needed when only eigenvalue approximations are desired.
+V = z0/norm(z0); 
+% Vectors storing entries of tridiagonal matrix \eqref{eq:tdcg}
+alph=zeros(k,1); bet = zeros(k,1);
+for j=1:k
+  q = A*V(:,j); alph(j) = dot(q,V(:,j));
+  w = q - alph(j)*V(:,j);
+  if (j > 1), w = w - bet(j-1)*V(:,j-1); end;
+  bet(j) = norm(w); V = [V,w/bet(j)];
+end
+bet = bet(1:end-1);

LectureCodes/Evp/lanczos/Matlab/lanczosev.m

+function [V,res] = lanczosev(A,k,j,w)
+V = [w/norm(w)];
+alpha = zeros(k,1);
+beta = zeros(k,1);
+res = [];
+for l=1:k
+  vt = A*V(:,l);
+  if (l>1)
+    vt = vt - beta(l-1)*V(:,l-1);
+  end
+  alpha(l) = dot(V(:,l),vt);
+  vt = vt - alpha(l)*V(:,l);
+  beta(l) = norm(vt);
+  V = [V, vt/beta(l)];
+  T2 = spdiags([[beta(1:l-1);0],alpha(1:l),[0;beta(1:l-1)]],[-1,0,1],l,l);
+  ev = sort(eig(full(T2)));
+  if (j > l)
+    ev = [zeros(j-l,1); ev]; 
+    res = [res; l ev']; 
+  else 
+    res = [res; l ev(end-j+1:end)'];
+  end
+end

LectureCodes/Evp/lanczos/Matlab/lanzevd.m

+% ex: lanzcosev
+n=100; M=gallery('tridiag',-0.5*ones(n-1,1),2*ones(n,1),-1.5*ones(n-1,1));
+[Q,R]=qr(M); A=Q'*diag(1:n)*Q; % eigenvalues \Blue{$1,2,3,\ldots,100$}
+ev = sort(eig(A));
+
+[V,res] = lanczosev(A,30,5,ones(100,1));
+figure; semilogy(1:30,abs(res(1:30,end-0)-ev(end-0)),'r-+',...
+		 2:30,abs(res(2:30,end-1)-ev(end-1)),'m-^',...
+		 3:30,abs(res(3:30,end-2)-ev(end-2)),'b-o',...
+		 4:30,abs(res(4:30,end-3)-ev(end-3)),'k-*');
+xlabel('{\bf step of Lanzcos process}','Fontsize',14);
+ylabel('{\bf |Ritz value-eigenvalue|}','fontsize',14);
+legend('\lambda_n','\lambda_{n-1}','\lambda_{n-2}','\lambda_{n-3}',3);
+
+print -depsc2 '../PICTURES/lanczosev1.eps';
+clear all;
+
+A = gallery('minij',100); ev = sort(eig(A));
+[V,res] = lanczosev(A,15,5,ones(100,1));
+figure; semilogy(1:15,abs(res(1:15,end-0)-ev(end-0)),'r-+',...
+		 2:15,abs(res(2:15,end-1)-ev(end-1)),'m-^',...
+		 3:15,abs(res(3:15,end-2)-ev(end-2)),'b-o',...
+		 4:15,abs(res(4:15,end-3)-ev(end-3)),'k-*');
+xlabel('{\bf step of Lanzcos process}','Fontsize',14);
+ylabel('{\bf error in Ritz values}','fontsize',14);
+legend('\lambda_n','\lambda_{n-1}','\lambda_{n-2}','\lambda_{n-3}',3);
+
+print -depsc2 '../PICTURES/lanczosev.eps';

LectureCodes/Evp/utils/eigensolversort.hpp

+#pragma once
 #include <Eigen/Dense>
 #include <vector>