基于LDL_T分解求实对称矩阵特征值的递归算法

２００８，４４（３）ＣｏｍｐｕｔｅｒＥｎｇｉｎｅｅｒｉｎｇａｎｄＡｐｐｌｉｃａｔｉｏｎｓ计算机工程与应用

１引言

矩阵特征值计算在许多工程领域中有着重要应用，因此，

多年来一直是个十分活跃的研究课题，并已研究出多种计算矩

阵特征值的算法［１－５］。这些算法各有优缺点，分别适用于不同场

合。本文基于线性代数和矩阵理论，将实对称矩阵的ＬＤＬＴ分

解用于估计该矩阵在给定区间内特征值的个数，进而拓展为一

种计算实对称矩阵特征值的递归算法。下面给出算法的推导与

理论分析，以及利用Ｍａｔｌａｂ对算法进行验证的结果，并对算法

的性能进行评价与测试。

２实对称矩阵的ＬＤＬＴ分解

由ＬＤＬＴ分解定理可知：若Ａ∈Ｒｎ×ｎ，ＡＴ＝Ａ，且Ａ的所有顺

序主子式都不为零，则存在唯一的单位下三角矩阵Ｌ和对角

矩阵Ｄ，使Ａ＝ＬＤＬＴ［６］。

为得到计算Ｌ、Ｄ元素的公式，设Ａ＝（ａｉｊ）ｎ×ｎ（其中ａｉｊ＝ａｊｉ），Ｌ

和Ｄ具有如下形式：

Ｌ＝１０…０

ｌ２１１…０

………

ｌｎ１ｌｎ２…

"

#

#

#

#

#

#

#

#

##

$

%

&

&

&

&

&

&

&

&

&&

’

１

，Ｄ＝

ｄ１１０…０

０ｄ２２…０

………

００…ｄｎｎ

"

(

(

(

(

(

(

(

(

((

$

%

&

&

&

&

&

&

&

&

&&

’

将上式代入Ａ＝ＬＤＬＴ，比较对应元素可得：

ａｉｊ＝

ｊ

ｋ＝１

)ｌｉｋｌｊｋｄｊｋ，ｊ＝１，２，…，ｎ；ｉ＝ｊ，ｊ＋１，…，ｎ

根据上式就可得到由ｊ＝１到ｊ＝ｎ逐列计算Ｌ、Ｄ各元素的公式为：

ｄ１１＝ａ１１

ｌｉ１＝ａｉ１

ａ１１

ｌｉｊ＝

ａｉｊ－

ｊ－１

ｋ＝１

)ｌｉｋｌｊｋｄｋｋ

ｄｊｊ

ｄｊｊ＝ａｊｊ－

ｊ－１

ｋ＝１

)ｌ２ｊｋｄｋｋ，ｊ＝１，２，…，ｎ；ｉ＝ｊ＋１，ｊ＋２，…，

*

,

,

,

,

,

,

,

,

,

,

,

+

,

,

,

,

,

,

,

,

,

,

,

-

ｎ

可见，上述算法保留了ＣｈｏｌｅｓｋｙＬＬＴ分解［６］的优点，而省去了开方运算。

３判断实对称矩阵在给定区间内特征值个数的算法为判断Ａ在区间（ａ，ｂ）内特征值的个数，根据实对称矩阵对角化原理［７］，特构造如下对称矩阵：

Ｐ（ａ）＝Ａ－ａＩ＝ＱＴｄｉａｇ（!１，!２，…，!ｎ）Ｑ－ａＱＴＱ＝ＱＴｄｉａｇ（!１－ａ，!２－ａ，…，!ｎ－ａ）Ｑ

Ｐ（ｂ）＝Ａ－ｂＩ＝ＱＴｄｉａｇ（!１，!２，…，!ｎ）Ｑ－ｂＱＴＱ＝ＱＴｄｉａｇ（!１－ｂ，!２－ｂ，…，!ｎ－ｂ）Ｑ

其中Ｑ∈Ｒｎ×ｎ，为正交矩阵；!１，!２，…，!ｎ为Ａ的特征值，且均为实数。

设矩阵Ｐ（ａ）和Ｐ（ｂ）的正惯性指数分别为Ｎ

ａ

＋

和Ｎ

ｂ

＋

，则由

惯性指数的定义可知［６］，Ｎ

ａ

＋

表示!１，!２，…，!ｎ中大于ａ的个数，

Ｎ

ｂ

＋

表示!１，!２，…，!ｎ中大于ｂ的个数，而Ｎ

ａ

＋

－Ｎ

ｂ

＋

即为矩阵Ａ在

基于ＬＤＬＴ分解求实对称矩阵特征值的递归算法

张鹍１，张有志２

ＺＨＡＮＧＫｕｎ１，ＺＨＡＮＧＹｏｕ－ｚｈｉ２

１．山东大学威海分校，山东威海２６４２０９

２．山东大学信息科学与工程学院，济南２５０１００

１．ＳｈａｎｄｏｎｇＵｎｉｖｅｒｓｉｔｙａｔＷｅｉｈａｉ，Ｗｅｉｈａｉ，Ｓｈａｎｄｏｎｇ２６４２０９，Ｃｈｉｎａ

２．ＳｃｈｏｏｌｏｆＩｎｆ．Ｓｃｉ．ａｎｄＥｎｇ．，ＳｈａｎｄｏｎｇＵｎｉｖｅｒｓｉｔｙ，Ｊｉ’ｎａｎ２５０１００，Ｃｈｉｎａ

Ｅ－ｍａｉｌ：ｚｈａｎｇｋｕｎｃｏｍ１１９＠ｓｉｎａ．ｃｏｍ

ＺＨＡＮＧＫｕｎ，ＺＨＡＮＧＹｏｕ－ｚｈｉ．ＲｅｃｕｒｓｉｖｅａｌｇｏｒｉｔｈｍｆｏｒｃａｌｃｕｌａｔｉｎｇｅｉｇｅｎｖａｌｕｅｓｏｆｒｅａｌｓｙｍｍｅｔｒｉｃｍａｔｒｉｘｂａｓｅｄｏｎＬＤＬＴｄｅｃｏｍｐｏｓｉｔｉｏｎ．ＣｏｍｐｕｔｅｒＥｎｇｉｎｅｅｒｉｎｇａｎｄＡｐｐｌｉｃａｔｉｏｎｓ，２００８，４４（３）：７８－８０．

Ａｂｓｔｒａｃｔ：ＡｒｅｃｕｒｓｉｖｅａｌｇｏｒｉｔｈｍｆｏｒｃａｌｃｕｌａｔｉｎｇｔｈｅｅｉｇｅｎｖａｌｕｅｓｏｆａｒｅａｌｓｙｍｍｅｔｒｉｃｍａｔｒｉｘｂａｓｅｄｏｎＬＤＬＴｄｅｃｏｍｐｏｓｉｔｉｏｎｉｓｇｉｖ－ｅｎ．Ｗｉｔｈｔｈｉｓａｌｇｏｒｉｔｈｍ，ｔｈｅｎｕｍｂｅｒｏｆｅｉｇｅｎｖａｌｕｅｓｏｆａｒｅａｌｓｙｍｍｅｔｒｉｃｍａｔｒｉｘｉｎｔｈｅｇｉｖｅｎｉｎｔｅｒｖａｌｃａｎｂｅｃｏｕｎｔｅｄ，ａｎｄｔｈｅｅｉｇｅｎ－ｖａｌｕｅｓｏｆｔｈｅｍａｔｒｉｘｃａｎｂｅｃａｌｃｕｌａｔｅｄ．Ｔｈｒｏｕｇｈｔｈｅｏｒｅｔｉｃａｌａｎａｌｙｓｉｓａｎｄｎｕｍｅｒｉｃａｌｓｉｍｕｌａｔｉｏｎ，ｔｈｉｓａｌｇｏｒｉｔｈｍｉｓｐｒｏｖｅｄｔｏｂｅｅｆ－ｆｅｃｔｉｖｅ．

Ｋｅｙｗｏｒｄｓ：ＬＤＬＴｄｅｃｏｍｐｏｓｉｔｉｏｎ；ｒｅａｌｓｙｍｍｅｔｒｉｃｍａｔｒｉｘ；ｅｉｇｅｎｖａｌｕｅ；ｒｅｃｕｒｓｉｖｅａｌｇｏｒｉｔｈｍ

摘要：基于线性代数与矩阵理论，给出利用ＬＤＬＴ分解计算实对称矩阵特征值的递归算法。该算法可求出实对称矩阵在给定区间内的特征值的个数，并可计算满足精度要求的特征值。理论分析和实际测试证明该算法是有效的。

关键词：ＬＤＬＴ分解；实对称矩阵；特征值；递归算法

文章编号：１００２－８３３１（２００８）０３－００７８－０３文献标识码：Ａ中图分类号：ＴＰ３０１．６；Ｏ１５１．２１

７８

２００８，４４（３）ｉ

０１２３４５６７ａｉ０．００００００．００００００．２５００００．２５００００．３１２５００．３４３７５０．３５９３８０．３５９３８ｂｉ

１．００００００．５０００００．５０００００．３７５０００．３７５０００．３７５０００．３７５０００．３６７１９ｂｉ－ａｉ

１．００００００．５０００００．２５００００．１２５０００．０６２５００．０３１２５０．０１５６２０．００７８１

ｉ

８９１０１１１２１３１４ａｉ

０．３５９３８０．３６１３３０．３６１３３０．３６１３３０．３６１３３０．３６１４５０．３６１５１ｂｉ

０．３６３２８０．３６３２８０．３６２３００．３６１８２０．３６１５７０．３６１５７０．３６１５７ｂｉ－ａｉ

０．００３９００．００１９５０．０００９７０．０００４９０．０００２４０．０００１２０．００００６表１

特征值计算过程

区间（ａ，ｂ）内的特征值的个数（含特征值的重数）。于是判断矩阵Ａ在区间（ａ，ｂ）内特征值个数的问题就转化为计算矩阵Ｐ（ａ）和Ｐ（ｂ）的惯性指数的问题。算法如下：

（１）对实对称阵Ｐ（ａ）进行ＬＤＬＴ分解得

Ｐ（ａ）＝Ｌ（ａ）Ｄ（ａ）ＬＴ（ａ）其中Ｌ（ａ）和Ｄ（ａ）分别为单位下三角阵和对角阵。显然，Ｌ（ａ）是非奇异的（其行列式为１）。

（２）计算对角阵Ｄ（ａ）的正惯性指数，即主对角线上大于０的元素个数；根据文献［６］中的定理１．１．３可知，Ｐ（ａ）和Ｄ（ａ）具

有相同的惯性指数，由此可求得Ｐ（ａ）的正惯性指数Ｎａ＋

。

（３）对Ｐ（ｂ）进行ＬＤＬＴ分解，求得Ｐ（ｂ）的正惯性指数Ｎｂ＋

；Ｎａ＋

－Ｎｂ＋

即为矩阵Ａ在区间

（ａ，ｂ）内的特征值个数。值得注意的是，当ａ或ｂ等于Ａ的某个特征值时，对Ｐ（ａ）

或Ｐ（ｂ）的ＬＤＬＴ分解将会失败，因为此时Ｐ（ａ）或Ｐ（ｂ）的顺序主子式会出现０。这时应适当调整原来指定的区间重新计算。

４计算实对称矩阵特征值的算法

设矩阵Ａ在指定区间（ａ，ｂ）内的特征值个数为ｋ，取区间

中点ｃ＝ａ＋ｂ２

，经计算得（ａ，ｃ）内的特征值个数为ｋ１，（ｃ，ｂ）内的

特征值个数为ｋ２，ｋ＝ｋ１＋ｋ２。若ｋ１＝０，则放弃区间（ａ，ｃ）；若ｋ１＞０则令ａ１＝ａ，ｂ１＝ｃ，再计算区间（ａ１，ｂ１）内的特征值个数。重复上述过程，直至区间（ａｉ，ｂｉ）内的特征值个数ｋｉ＞０且｜ｂｉ－ａｉ｜＜ε（其中ε为预定计算精度），则取ａｉ＋ｂｉ２

作为矩阵Ａ的一个特征值，其重

数近似为区间（ａｉ，ｂｉ）内的特征值个数ｋｉ。上述算法的流程如图

１所示。

要计算区间（ａ，ｂ）内的所有特征值，就要遍历对分得到的

所有包含特征值的子区间。记录这些子区间的端点是极其繁琐的。为此，可利用函数递归调用的方式自动保存这些参数。设利用ＬＤＬＴ分解求对称矩阵Ｐ（ｘ）＝Ａ－ｘＩ的正惯性指数的函数为ｎｘ＝ｉｎｅｒｔｉａ＿ｐｏｓ（Ａ，ｘ），则计算实对称矩阵Ａ在区间（ａ，ｂ）内所有特征值的递归算法如下：

ｆｕｎｃｔｉｏｎｅｉｇ＿ｌｄｌ

（Ａ，ａ，ｂ，ｎａ，ｎｂ，ε）ｃ＝（ａ＋ｂ）／２／／区间对分

ｉｆ｜ｂ－ａ｜＜εｔｈｅｎ／／精度达到要求

ｆｏｒｉ＝１ｔｏ｜ｎａ－ｎｂ｜

ｏｕｔｐｕｔｃ／／输出特征值，重数为｜ｎａ－ｎｂ｜ｅｌｓｅ／／精度未达到要求

ｎｃ＝ｉｎｅｒｔｉａ＿ｐｏｓ

（Ａ，ｃ）／／计算区间中点的正惯性指数ｉｆｎａ≠ｎｃｔｈｅｎ／／

（ａ，ｃ）内有特征值ｃａｌｌｅｉｇ＿ｌｄｌ

（Ａ，ａ，ｃ，ｎａ，ｎｃ，ε）／／递归，计算（ａ，ｃ）内的特征值ｉｆｎｂ≠ｎｃｔｈｅｎ／／

（ｃ，ｂ）内有特征值ｃａｌｌｅｉｇ＿ｌｄｌ

（Ａ，ｃ，ｂ，ｎｃ，ｎｂ，ε）／／递归，计算（ｃ，ｂ）内的特征值ｒｅｔｕｒｎ

在函数ｅｉｇ＿ｌｄｌ中将区间端点的正惯性指数作为参数记忆

和传递，是为了避免递归时在同一端点上重复进行ＬＤＬＴ分解。上述算法保证每次调用产生新的、必要的子区间，并且在新产生的两个子区间的公共端点处只进行一次ＬＤＬＴ分解。

主程序中需要预先计算矩阵Ｐ（ｘ）在初始区间端点处的正惯性指数，算法如下：

ｆｕｎｃｔｉｏｎｅｉｇ＿ｌｄｌ＿ｍａｉｎ（Ａ，ａ，ｂ，ε）ｎａ＝ｉｎｅｒｔｉａ＿ｐｏｓ（Ａ，ａ）ｎｂ＝ｉｎｅｒｔｉａ＿ｐｏｓ

（Ａ，ｂ）ｃａｌｌｅｉｇ＿ｌｄｌ（Ａ，ａ，ｂ，ｎａ，ｎｂ，ε）ｒｅｔｕｒｎ

可见，实现上述递归算法的程序是非常简洁的。根据矩阵特征值估计理论［８］，利用上述算法可以求得实对称矩阵Ａ的所

有特征值。

５算法性能分析及测试验证

５．１算法的复杂度

设Ａ∈Ｒｎ×ｎ

，初始区间为（ａ，ｂ），精度为ε。采用对分法找到

（ａ，ｂ）内的某个特征值理论上需要对初始区间分割Ｎ＝［ｌｂｂ－ａ

ε

］＋

１次，进行ＬＤＬＴ分解的次数约为［ｌｂｂ－ａ

ε

］＋３次

（包括端点ａ，ｂ）。

可见，精度每提高一个数量级，只需增加３￣４次对分。例如，随机取对称矩阵：

Ａ＝４２０２２１０５９０５５４２９４#

$$$$$$$$$$%

&’

’

’’

’

’

’

’’

’(

９其特征值为０．３６１５４，１．９４４０９，４．２４５１６，２１．４４９２０。

若用上述算法在（０，１）内计算特征值，精度取１０－４。计算过

程如表１所示（其中Ｎａ＋－Ｎｂ＋

＝１）。

由表１可知，最终得到的特征值为（ａ１４＋ｂ１４）／２＝０．３６１５４，与实际值相比达到指定精度。理论上的分割次数为Ｎ＝［ｌｂ１０４］＋１＝

１４，实际分割次数也是１４。

若要求出（ａ，ｂ）内的所有ｋ个特征值，最坏情况下约需ｋＮ

给ａ，ｂ，ε赋初值对Ｐ（ａ），Ｐ（ｂ）做ＬＤＬＴ分

解求得Ｎａ＋

和Ｎｂ

＋

调整ａ，ｂ

分解失败？

Ｎａ＋

＝Ｎｂ＋

？

｜ｂ－ａ｜＜ε

？令ａ＝（ｂ＋ａ）／２或ｂ＝（ａ＋ｂ）／２

特征值"≈（ａ＋ｂ）／２重数为Ｎａ＋

－Ｎｂ

＋

Ｙ

Ｙ

ＮＮ

Ｎ

Ｙ

图１算法流程图注：调整ａ，ｂ的方法是：若Ｐ（ａ）分解失败则令ａ＝ａ－δ（δ充分小）；若Ｐ（ｂ）分解失败则令ｂ＝ｂ＋δ。若Ｎａ＋

＝Ｎｂ＋

（即该区间内无特征值）则更换为由对分生成的其它子区间。

张鹍，张有志：基于ＬＤＬＴ分解求实对称矩阵特征值的递归算法

７９

２００８，４４（３）ＣｏｍｐｕｔｅｒＥｎｇｉｎｅｅｒｉｎｇａｎｄＡｐｐｌｉｃａｔｉｏｎｓ计算机工程与应用

双音节词这些产业类型母亲法语四万只有成为显示衰老平均值

传统算法３．５３．０３．０４．５２．５２．５３．５２．５４．５４．０３．３５

本文算法３．５４．０４．０３．５３．５２．５３．５４．０４．５４．０３．７

传统算法１．９９

－１．５６３．９８－３．２６－３．９４－１．３１０．３０３．６４０．３１－１．４６－０．１３１

本文算法

３．７６１．４０４．１９－５．６２６．８６２．３１１．３７３．６５０．４０－０．４９１．７８３

ＭＯＳＳＮＲ／ｄＢ

表１两种频谱平滑方法的ＭＯＳ得分和信噪比（ＳＮＲ）

以看到，共振峰的不连续情况得到了较大的改善，音节间的共振峰实现了较为平滑的过渡。

为了与传统的频谱平滑算法进行比较，也使用了基于语音的源－滤波器模型的频谱平滑算法对这１０个双音节词进行了实验。并将这２种方法在ＭＯＳ得分和信噪比（ＳＮＲ）方面进行了对比，如表１所示。

４结论

通过在主观听音和客观评价两方面与传统算法的比较，频

谱平滑算法显示出了较好的特性，使合成语音的频谱得到了较好的修正，同时其音质没有明显的下降。

然而，由于算法需要进行ＦＦＴ计算，因此，其计算量比传统方法有所增加。

在语音频谱平滑算法的研究中还有很多问题需要更好的解决方法。例如：如何找到一种更加准确而简洁的方法对语音的频谱进行表示，如何更加自然的从源语音帧的频谱向目标语音帧的频谱进行过渡等，这些问题是今后进一步研究的方向。（收稿日期：２００７年８月）

参考文献：

［１］ＱｕａｔｉｅｒｉＴＦ．Ｄｉｓｃｒｅｔｅ－ｔｉｍｅｓｐｅｅｃｈｓｉｇｎａｌｐｒｏｃｅｓｓｉｎｇ：ｐｒｉｎｃｉｐｌｅｓａｎｄ

ｐｒａｃｔｉｃｅ［Ｍ］．Ｂｅｉｊｉｎｇ：ＰｕｂｌｉｓｈｉｎｇＨｏｕｓｅｏｆＥｌｅｃｔｒｏｎｉｃｓＩｎｄｕｓｔｒｙ，２００４．［２］ＰｆｉｔｚｉｎｇｅｒＨＲ．ＤＦＷ－ｂａｓｅｄｓｐｅｃｔｒａｌｓｍｏｏｔｈｉｎｇｆｏｒｃｏｎｃａｔｅｎａｔｉｖｅ

ｓｐｅｅｃｈｓｙｎｔｈｅｓｉｓ［Ｃ］／／ＰｒｏｃＩＣＳＬＰ２００４，Ｋｏｒｅａ，２００４：１３９７－１４００．

［３］ＷｏｕｔｅｒｓＪ，ＭａｃｏｎＭＷ．Ｃｏｎｔｒｏｌｏｆｓｐｅｃｔｒａｌｄｙｎａｍｉｃｓｉｎｃｏｎｃａｔｅｎａ－

ｔｉｖｅｓｐｅｅｃｈｓｙｎｔｈｅｓｉｓ［Ｊ］．ＩＥＥＥＴｒａｎｓｏｎＳｐｅｅｃｈａｎｄＡｕｄｉｏＰｒｏｃｅｓｓ－

ｉｎｇ，２００１，９

（１）：３０－３８．［４］ＫａｎｇＨ，ＬｉｕＷｅｎｊｕ．Ｓｉｎｕｓｏｉｄａｌａｌｌ－ｐｏｌｅｍｏｄｉｆｉｃａｔｉｏｎｂａｓｅｄｓｐｅｃｔｒａｌ

ｓｍｏｏｔｈｉｎｇｆｏｒｃｏｎｃａｔｅｎａｔｉｖｅｓｐｅｅｃｈｓｙｎｔｈｅｓｉｓ［Ｃ］／／ＮａｔｕｒａｌＬａｎｇｕａｇｅ

ＰｒｏｃｅｓｓｉｎｇａｎｄＫｎｏｗｌｅｄｇｅＥｎｇｉｎｅｅｒｉｎｇ，ＩＥＥＥＮＬＰ－ＫＥ’

０５，２００５．［５］ＷｏｕｔｅｒｓＪ，ＭａｃｏｎＭＷ．Ｓｐｅｃｔｒａｌｍｏｄｉｆｉｃａｔｉｏｎｆｏｒｃｏｎｃａｔｅｎａｔｉｖｅ

ｓｐｅｅｃｈｓｙｎｔｈｅｓｉｓ［Ｃ］／／Ａｃｏｕｓｔｉｃｓ，Ｓｐｅｅｃｈ，ａｎｄＳｉｇｎａｌＰｒｏｃｅｓｓｉｎｇ２０００ＩＣＡＳＳＰ’００，Ｐｒｏｃｅｅｄｉｎｇｓ２０００ＩＥＥＥＩｎｔｅｒｎａｔｉｏｎａｌＣｏｎｆｅｒｅｎｃｅ，２０００，２：９４１－９４４．

［６］蔡莲红，黄德智，蔡锐．现代语音技术基础与应用［Ｍ］．北京：清华大

学出版社，２００３．

［７］清华大学计算机科学与技术系．语音合成语料库ＴＨ－ＣｏＳＳ技术报

告［Ｒ］．北京：清华大学，２００３．

（上接７１页）

次分割（与特征值在该区间内的分布有关，特征值越密集或有重特征值时每次迭代划分的子区间越少）。可以证明，当矩阵阶数ｎ足够大时，一次ＬＤＬＴ分解约需ｎ３／６次乘除法。每次分割要进行一次ＬＤＬＴ分解，因此，求出所有ｋ个特征值需要的乘除法次数大约为ｋＮｎ３／６。

５．２算法的有效性从算法的推导可以看出，该算法可以完成指定精度的特征值的计算。在对称阵Ａ存在重特征值或特征值比较密的情况下，为得到较准确的结果，应该保证足够高的精度。从上面的分析可知，算法的复杂度随精度的提高呈对数特性，故提高精度是可行的。

例如，求Ａ＝ＱＴｄｉａｇ（－１，０，１．０００１，１，１，２）Ｑ（Ｑ为一正交阵）的所有特征值。取初始区间为（－２，３），精度为１０－３，利用上述算法求得特征值为－０．９９９９，－０．０００２，１．０００１，１．０００１，１．０００１，

１．９９９９。在此情况下，算法判断Ａ有３重特征值１．０００；若精度提高至１０－４，则得特征值为－０．９９９９８，０．００００１，０．９９９９９，０．９９９９９，

１．００００７，１．９９９９８。

此时算法判断Ａ有２重特征值１．００００和１重特征值１．０００１。

可见，只要预定精度足够高，该算法就可以求得满足精确要求的特征值。究竟取多高的精度才算合适，应根据实际需要来决定。

５．３算法的稳定性

该算法中，特征值的计算只用到了加法和除以２的运算，运算简单，对误差不敏感；ＬＤＬＴ分解虽然计算相对复杂，但由于只需对对角阵Ｄ中对角线元素的符号进行判断，分解过程中的误差对特征值的计算影响不大。

当区间端点靠近实际特征值时，ＬＤＬＴ分解得到的对角阵

Ｄ的某些对角元素绝对值很小，在有限精度的情况下容易判断为０，造成惯性指数的计算错误。这时应略微扩展原区间重新计算。

例如，计算Ａ＝ＱＴｄｉａｇ（－１，０，１，１，２）Ｑ（Ｑ为正交矩阵）在（０，２）内的特征值时ＬＤＬＴ分解失败。此时只需将原区间扩展

为（１０－６，２＋１０－６），即可得到正确结果。

理论和实验表明，若成功完成ＬＤＬＴ分解，就能得到满足给定精度要求的结果；当ＬＤＬＴ分解失败时，微调区间端点，也能获得正确结果。

６结论

本文设计的利用ＬＤＬＴ分解计算实对称矩阵特征值的递归算法能够判断实对称矩阵在给定区间内特征值的个数，并能计算该区间内的特征值。算法实现简单，稳定性好，提高精度需要的附加运算也不多。（收稿日期：２００７年７月）

参考文献：

［１］魏立峰，李晓梅．计算实对称矩阵广义特征值问题的并行算法［Ｊ］．计

算机工程与应用，２００１，３７（１１）：４－５．

［２］ＨａｓａｎＭＡ．Ｎｅｗａｌｇｏｒｉｔｈｍｓｆｏｒｃｏｍｐｕｔｉｎｇｔｈｅｍｉｎｉｍｕｍｅｉｇｅｎｐａｉｒ

ｏｆｔｈｅｇｅｎｅｒａｌｉｚｅｄｓｙｍｍｅｔｒｉｃｅｉｇｅｎｖａｌｕｅｐｒｏｂｌｅｍ［Ｃ］／／２００２ＩＥＥＥＩｎ－ｔｅｒｎａｔｉｏｎａｌＳｙｍｐｏｓｉｕｍｏｎＣｉｒｃｕｉｔｓａｎｄＳｙｓｔｅｍｓＰｒｏｃｅｅｄｉｎｇｓ，ＩＥＥＥＰｒｅｓｓ，Ｐｈｏｅｎｉｘ，ＵＳＡ，２００２，４：７６７－７７０．

［３］ＨｅＤａ－ｋｅ，ＷａｎｇＪｉａｎ－ｂｏ．Ｐａｒａｌｌｅｌｃｏｍｐｕｔｉｎｇｏｆｅｉｇｅｎｖａｌｕｅｏｆｄｏｕ－

ｂｌｙｓｔｏｃｈａｓｔｉｃｍａｔｒｉｘ［Ｃ］／／ＰｒｏｃｅｅｄｉｎｇｓｏｆＦｉｆｔｈＩｎｔｅｒｎａｔｉｏｎａｌＣｏｎｆｅｒ－ｅｎｃｅｏｎＡｌｇｏｒｉｔｈｍｓａｎｄＡｒｃｈｉｔｅｃｔｕｒｅｓｆｏｒＰａｒａｌｌｅｌＰｒｏｃｅｓｓｉｎｇ．Ｂｅｉ－ｊｉｎｇ，Ｃｈｉｎａ：ＩＥＥＥＰｒｅｓｓ，２００２：３５５－３５８．

［４］ＬｉｕＹｉ－ｇｕａｎｇ，ＹｏｕＺｈｉ－ｓｈｅｎｇ，ＣａｏＬｉ－ｐｉｎｇ．Ａｃｏｎｃｉｓｅｆｕｎｃｔｉｏｎａｌ

ｎｅｕｒａｌｎｅｔｗｏｒｋｃｏｍｐｕｔｉｎｇｔｈｅｌａｒｇｅｓｔ

（ｓｍａｌｌｅｓｔ）ｅｉｇｅｎｖａｌｕｅａｎｄｏｎｅｃｏｒｒｅｓｐｏｎｄｉｎｇｅｉｇｅｎｖｅｃｔｏｒｏｆａｒｅａｌｓｙｍｍｅｔｒｉｃｍａｔｒｉｘ［Ｃ］／／Ｐｒｏｃｅｅｄ－ｉｎｇｓｏｆ２００５ＩｎｔｅｒｎａｔｉｏｎａｌＣｏｎｆｅｒｅｎｃｅｏｎＮｅｕｒａｌＮｅｔｗｏｒｋｓａｎｄＢｒａｉｎ，ＩＥＥＥＰｒｅｓｓ，Ｂｅｉｊｉｎｇ，Ｃｈｉｎａ，２００５：１３３４－１３３９．

［５］ＹａｎｇＬｉｕ，ＢｏｕｇａｎｉｓＣＳ，ＣｈｅｕｎｇＰＹＫ，ｅｔａｌ．Ｈａｒｄｗａｒｅｅｆｆｉｃｉｅｎｔ

ａｒｃｈｉｔｅｃｔｕｒｅｓｆｏｒｅｉｇｅｎｖａｌｕｅｃｏｍｐｕｔａｔｉｏｎ［Ｃ］／／ＰｒｏｃｅｅｄｉｎｇｓｏｆＤＡＴＥ’０６．Ｍｕｎｉｃｈ，Ｇｅｒｍａｎｙ：ＳｉＧＰｕｂｌｉｓｈ，２００６，１：１－６６．

［６］蔡大用，白峰杉．现代科学计算［Ｍ］．北京：科学出版社，２０００．［７］熊全淹，叶明训．线性代数［Ｍ］．３版．北京：高等教育出版社，１９８７．［８］程云鹏．矩阵论［Ｍ］．西安：西北工业大学出版社，２０００．

８０