运用递归的方法求解对称三对角矩阵的特征值
问题
求解如下的对称三对角矩阵的特征值



我们在求解过程中，发现有如下的递推关系



于是我想到了运用matlab中的sym变量并用递归的方法把特征多项式表达出来，然后求解。
递归函数
首先我们定义这样一个函数，理解不畅可以参考斐波那契数列的递归求解
function y = recurMatrix( n )
%recurMatrix 运用递归的方法求解对称三对角矩阵特征值
%   n：size of the matrix,n>0 and is int
%   y：特征多项式
syms y lambda;
if n == 1
    y = lambda-2;
elseif n==2
    y = (lambda-2)^2-2;
else 
    y = (lambda-2)*recurMatrix(n-1) - recurMatrix(n-2);
end
然后我们可以输入任意的n，运用solve函数进行求解
效果

结果出来的是分数表达式如果想要具体的数值，可以使用eval/double


写在最后
当然了，我所解决的问题依旧是矩阵规模较小的情况，而且矩阵具有比较强的特殊性。
这个问题可以当成训练递归思想的一个例子，但在实际问题中，如果真的求解一般矩阵的特征值，那还得是诸如QR方法之类的比较常见的，成熟的方法

前言：以下内容不是严格的数学表述， 以自己理解的思路形式叙述。
二次型：
这个名词是来自于线性代数， 多用于二次规划和优化组合等问题。
在线性代数里形如以下函数表达式称为二次型：（A是对称矩阵）
这里Q(x)输出的是一个常数。
(输入不同的x时，输出不同的常数，在n+1维空间中展示图形，如果x是二维的向量， 可以以三维来表示所有二维向量经过矩阵A变换后的样子）
下图是 
 中的向量x，对应z的图形展示， 变成三维图象：
根据个人理解，只要是方程等式右边理解为变量时， 在图形上就要多出一个维度才可以表示。所以，以上表达式可以在n+1维空间中展示。
根据对角化分解， 
 , 其中，P是特征向量组成的矩阵， D是对角线上为特征值的矩阵。
根据对称矩阵分解定理， 当 
 ，时 
 中的特征向量正交。
定义及应用：
一个二次型Q是：
a, 正定的，如果 
  , 有 
 .      （存在最小值）
b, 负定的，如果 
  , 有 
 .      （存在最大值）
c, 不定的， 如果 
 的值 ，即有正， 又有负。  
d, 半正定的，如果 
  , 有 
 .   （最小值为0）
e, 半负定的，如果 
  , 有 
 .    （最大值 为0）
分析和理解二次型：
根据二次型公式， 可以理解为矩阵A把向量的空间进行变换， 即对主轴进行了变换（有可以是缩放、旋转、翻转等）， 变换后的向量再映射成一个实数
对称矩阵的特点， 可以分解成 
 , 其中D是对角矩阵， 对角线上的是特征值：
有N个特征值 和特征向量
这种矩阵对应的是旋转和拉伸， 没有压缩， 所以可以分析对那一个特征向量进行了拉伸。
根据矩阵点积的意义，如果 A = 
 可以把公式分解成：
该公式正是 
, 模长的平方，是一个二次多项式。
为什么？因为该公式的代数表达式就是一个多元二次方程。
注意这里的表达式中，分析的目标是矩阵A， 一定要注意。
任意维度的空间都可以有二次型。
任意对称矩阵都可以有二次型形式.
二次型分解的几种情况：
一， 当二次型中的矩阵A是对称矩阵，但对角线不是特征值时， 需要进行变量代换， 如把A分解成 
 。
设:      
这时二次型   
对于 
 时, 最大值和最小值是特征值的最大值和最小值。 

关联知识

Python
Numerical

介绍
雅可比法
相似变换和对角化
雅各比旋转
雅可比对角化
转换为标准形式
幂和逆幂方法
逆幂方法
特征值漂移
幂方法
Householder简化为对角线形式
Householder矩阵
对称矩阵的Householder约简
累积转换矩阵
Sturm序列
格什戈林定理
括号化特征值
特征值的计算
特征向量的计算
对称三对角矩阵的特征值
…更多
详情参阅http://viadean.com/python_numerical_matrix.html

 Lanczos法的目的：将实对称矩阵转换成为三对角矩阵（稀疏矩阵），从而便于计算机储存和后续的计算。

在三对角矩阵矩阵上，采用QR分解，得到矩阵的特征值。

# -*- coding: utf-8 -*-
"""
Created on Tue Nov 20 12:41:40 2018

@author: yujin.wang
"""
import numpy as np
import scipy as sci
import matplotlib.pyplot as plt
import time
import sys

def lanczos(A,b,nmax):
	m = np.size(A,1)
	alpha = []
	beta = [0]
	qprev = np.zeros(m)
	q = b / np.linalg.norm(b)
	for n in range(nmax):
		v = np.dot(q,A) #*sci.linalg.inv(q)
		temp = np.dot(v,q.T)
		alpha.append(temp[0][0,0])
		v = v - np.dot(beta[-1],qprev)-np.dot(alpha[-1],q)
		beta.append(sci.linalg.norm(v))
		qprev = q
		q = v/beta[-1]
	beta = beta[1:-1]
	T = np.diag(alpha) + np.diag(beta,1) +np.diag(beta,-1)
	return T


def qreigen(A,num=100):
	m = np.size(A,1)
	p = np.eye(m)
	for i in range(num):
		v = np.diag(A)
		q,r = sci.linalg.qr(A)
		A = np.dot(r,q)
		p = np.dot(p,q)
		tol = max(np.diag(A))-max(v)
		print 'TOL:',tol,'Max. Eig:',max(v)
		if np.abs(tol) <1e-6:
			break
#	s = np.diag(np.diag(r))
	return p,np.diag(A)

if __name__ == '__main__':
	A = np.matrix([[5,1,3],[1,2,0],[3,0,11]])
	print np.linalg.eig(A)[0]
	print '$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
	p,s = qreigen(A)
	print s
	print '$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
	b = np.matrix([1,1,1])
	tridiag = lanczos(A,b,23)
	p,s = qreigen(tridiag,num=200)
	print max(s),min(s)

参考：https://baike.baidu.com/item/Lanczos%E7%AE%97%E6%B3%95/9849921?fr=aladdin 

           http://www.cnblogs.com/qxred/p/dcalgorithm.html

           http://www.doc88.com/p-1446190703774.html