D0I:10.13374/j.issn1001-053x.1987.01.008 北京钢·铁学院学报 第9卷第1期 Journal of Beijing University Val.9 No.1 1987年1月 of Iron and Steel Technology Jan。1987 AR模型的两种 计算量较小的建模方法 张思宇 陈克兴 (冶金机械教研室) 摘 要 建棋方法是时间序列分析的核心问题之一，本文给出了两种基子最小二乘法的自回归模 型(AR模型)的建棋方法，采取预留少量数据递补进入计算的办法，使矩阵XTX可以用分块 矩阵求逆公式递推求逆，或者用矩阵的Crouts分解法递推求解。同时引入了Wino唱rad向盘内 积快速算法，充分利用各向量和各矩阵之间的关系来减少计算工作量。使计算量比一般最小 二乘建棋方法大幅度减少，达到与Marple©算法和Butg的最大熵谱法可比的径度， 关键词：建模，数据递补，递推 Two Modelling Methods of AR with Less Computational Operations Zhang Siyu Chen Kexing Abstract Method of modelling is a major question of the time series analysis.Two modelling methods,of autoregressive (AR)model based on the LS approach are proposed in this paper.By means of reserving a few data to enter calcu- lating successively,the matrix XTX can be inversed recursively by the inver- sing formula of block matrix or result from the Crout resolving method.Mo- reover,a quick scalar product algorithm,advanced by S.Winograd,is emp- 1985-10-17收税 48
第 卷 第 期 北 京 钢 铁 学 院 学 报 入 。 年 月 。 模型的两种 计算量较小的建模方法 张思宇 陈克兴 冶金机械教研室 摘 要 建模方法是时 间序列分析的核心 问题之一 本文给出了两种基于最小二乘法的 自回归 模 型 人 模型 的建模方法 。 采取预留少量数据递补进入计算的办法 ， 使矩阵 可以 用分块 矩阵求逆公式递推求逆 ， 或者 用矩阵的 分解法递推求解 同时引入了 向量 内 积快速算法 ， 充分利用各向量和各矩阵之间 的关系来减少计算工作量 。 使计算量比一般最小 二乘建模方法大幅度减少 ， 达 到与 算法和 的最大嫡谱法可比的程度 。 关键词 建模 ， 数据递补 ， 递推 夕 少“ ” 夕 衫 价 ， ， ， ， 一 一 收 稿 DOI ：10．13374／j ．issn1001－053x．1987．01．008
loyed and the relations between vectors or matrices are applied for reducing computational operations,The operations have been largely reduced so that they are less than the operations of normal LS approach,and they can be compared with the famous Marple algorithm and Burg's maximum entropy spectral analysis method. Key words:modelling,supplement data successively,recursion 引 言 自回归模型(AR模型)是时间序列分析方法常采用的模型。它的形式简单，建模比较容 易。在许多场合还用来代替较为复杂的自回归滑动平均模型(ARMA模型)，以简化模型结 构和减少建模计算工作量1)。 建模方法是时间序列分析的核心问题之一。AR模型的建模方法有最小二乘估计、最小 平方和估计、极大似然估计等几类(23)。其巾基于最小二乘原理的建模方法有多种，并且包括 几种递推算法(3,4)。本文基于最小二乘估计的基本原理，在样本数据的取用上与常用的AR 建模方法不同，从而导出了两种新的建模方法。其中引入了Winograd向量内积快速算法， 并充分利用建模过程中各量之间的关系来减少计算工作量，使之达到与Bu「g的最大熵谱法和 Marple算法(5)可比的程度。 1预留候补数据和递推求逆 设有一组样本数据{X,X,…,X},它构成一个平稳的、均值为零的时间序列。可以 用一个n阶的自回归模型来描述它： X:=Pa,IXt-1+a,zXt-2+..+aXt-n+at a:~NID(0,o)。 (1) 这个模型记为AR(（n)。建模的任务是以一定的方法估计出上式中的n个模型系数(Pn,, P,2,…,P。,。)。最小二乘估计定义残差平方和 E,,6=之a=之(X-pa,1X1-po,eX-2--Pain Xt-n)2。（2） t-a 这里若改变，b可以得到不同的数据取用形式。然后，利用多元函数求极值的方法，在残差平 方和最小的意义下估计出AR模型的系数。用矩阵表示为Φ，=(XtX)-1XY。当a= n+1,b=N时(N是样本长度)，上式中 Xn………X X+1 X= Pa,n XN-1…XN-n XN 49
五 ， ， ， ， 引 言 自回归模型 模型 是时 间序列分析方法常采 用的模型 。 它 的形式简单 ， 建模比较容 易 。 在许多场 合还 用 来代替较为 复杂 的 自回归滑动平均模型 模型 ， 以简化模型 结 构 和减少 建模计 算工作量〔 〕 。 建模方法是 时 间序列分析的核心 问题之一 。 模型 的建模方法有最小二 乘估 计 、 最 小 平方和估计 、 极大似然估计 等几 类， 〕 。 其 中基 于最小二乘原理 的 建模方法有 多种 ，并 且包 括 几种递 推算法 〔 ” ， 魂〕 。 本文基于最小二 乘估计的基本原理 ， 在样本数据的取 用 上与常用 的 建模方法 不 同 ， 从而 导 出 了两 种新的建模方法 。 其 中引人 了 向量 内积快速算法 ， 并充分利 用 建模过程中各量之 间的关 系来减少计 算工作量 ， 使 之达到与 的最大嫡谱法 和 算法〔 〕可 比的程度 。 预 留候补数据和递推求逆 设 有一 组样本数 据笼 ， ， … ， ， 它构 成一个平稳 的 、 均值为 零的时 间序列 。 可 以 用 一个 阶的 自回归模型 来描述它 甲 。 ， 一 甲 。 ， 一 … 甲 。 ， 。 一 。 ， 孟 。 这 个模型记为 。 建模的任 务是 以一定 的方法 估计 出上 式 中的 个模型系数 甲 。 ， ，… ，甲 。 ， 。 。 最 小二乘估计 定 义残 差平方和 ‘ 甲 ， ， 。 ， ， 万 圣 艺 一 甲 。 ， ， 一 ， 一 甲 。 ， 卜 一 、一 甲 。 ， 。 卜。 。 这 里 若改 变 ， 可 以 得到不 同的数 据取用 形 式 。 然后 ， 利用 多元 函数 求极值的方法 ，在残差平 方 和最小 的意 义 下估计 出 模型 的系数 。 用 矩 阵 表 示 为巾 ‘ 二 一 ’ 丁 。 当 二 ， 时 是 样本长度 ， 上式中 示 二 一 · · … 一 。 一 、 刁‘ … ， … ︿切甲
这就是常用的普通最小二乘AR模型建模方法。 最小二乘估计具有精度高，所得到的自回归谱的性能优良等优点。如文献〔3〕指出的， 它是目前最基本的AR模型的建模方法之一。然而上述公式的直接求解，无论采用什么解法计 算量都比较大；并且它一般只适合一次去估计一个某阶数的模型，在最佳模型的阶数未知的 情况下是不利的。这是普通最小二乘建模方法的两个主要弱点。这里采用预留少量数据递补 进入计算的办法，使这两个弱点同时得到改善，并且不损失估计精度。 令样本数据序列的排列为{X-m,X-m+1,…,X一1，X0,X:,…,XN},其中{X-m, X-m+1,…X一1}一段称为候补数据序列。 令a=1,b=N,由(2)式对于AR(n)有： E,m=芝a,=芝(X,-9a,1-1--paaX-a) (3) 它的解Dn=Cpa,:,pn,2,…,pa,n)T=（XXa）-1XY, (4) 其中： Xo X,…X-n+1 X0…X-n+2 X: Xn= (5) XN-1 XN-2 XN- X 这里数据的取用长度是{X一n十1，·，X、}。每次升阶时，由候补数据序列中依次递补上 个数据。对于AR(n+1)模型有 E+1,1,N=空a:=芝（X,-pa+1,1X-1--pa+1,+1X-a-） (6) t=1 t=1 解为： n+1=[+1,,n+1,n1=(XTo+iXn+1)-1X'n+iY (7) 其中 X0X1…X-m+1|X-n X1 X:X0…X-n+2|X-a+1 X2 Xm+1= 」： =〔XnZn〕，Y= (8) l8 XN-1 XN-2 XN-n XN-n-1 N 这里数据的取用长度是{X-,…,XN}。递补了X-n一个数据。X矩阵的行数不变，右侧 增加了一列Z。,Y矩阵不变。这时有性质： XX+-〔2]x2-〔22z, (9) 4v-[)Y=CY. (10) 这里XTn+1Xa+1矩阵的左上角n×n阶方阵就等于X:Xm,且XTa十1Y列向量的前n个元素 就等于XY。(9)式的结构刚好可以利用矩阵分块求逆公式来计算它的逆，并且与低一阶 50
这就是常用 的普通最小二乘 模型建模方法 。 最小二乘估计 具有精 度高 ， 所得到 的 自回归谱的性能优 良等优 点 。 如文献 〔 〕 指出的 ， 它是 目前最基 本的 模型的建模方法 之一 。 然而上述 公式的直接求解 ， 无论采 用什么 解法计 算量都 比较大 并且它 一般 只适合一次去估计一 个某阶数 的模型 ， 在最佳模型 的阶数未知 的 情况 下是 不利 的 。 这 是 普 通 最小 二 乘建模方法 的两 个主要 弱 点 。 这里 采 用预 留少量数据递 补 进 人计算的办法 ， 使这两 个弱 点 同时得到 改善 ， 并且不损失估计精度 。 令样本数据 序列 的排列为 一 ， 一 ， … ， 一 ， 。 ， ， … ， ， 其中 一两 一 ， … 一 一 段称 为候补数 据序列 。 令 ， ， 由 式对 于 有 ， ， ，。 二 艺 艺 一 甲 。 ， 一 一 · 一 甲 ， 。 一 一 盛 它 的解 其 中 巾 。 〔切 。 ， ， 甲 ， ， … ， 印 。 ， 。 〕 百 。 一 百 ， … 一 … 一 悦 ，﹄ 一 … 。 一 这 里 数 据的取用长度是 一。 ， … ， 。 每 次升 阶时 ， 由候补数据序列 中依次递 补上 一个数据 。 对 于 模型有 ， ， ， 、 艺 圣 叉 一 尸、 ， 甲 一 甲 ， ， ， 一 · 一 甲 ， 。 一 一 “ 解 为 其 中 中 〔 印 ， ，， … 〕 二 了 一 叫 一 、 十 气 一 一’ 二 一 一 一 一 一 ， 〔 ‘ 〕 ’ ‘ 七…、 … ﹄ 这 里数据 的取 用长度 是 一 。 ， … ， 。 递 补了 石 一个数据 。 矩 阵的行数不 变 ， 右侧 增加 了一列 。 ， 矩 阵不 变 。 这 时有性质 ’ 一 十 〔和 〔 · · 〕 〔数 尝〕 〔 撰〕 〔 粼〕 。 这 里 叶 矩 阵的左 上 角 只 渝方 阵就等 于 工 。 ， 且 。 列 向量 的前 个元 素 就等于 万 。 式的结构刚好 可 以利 用矩阵分块求逆公式来计算它 的逆 ， 并且与低一阶
矩阵的逆(XXn)1构成递推关系。(10)式则直接给出了XT。+1Y与XY的关系。 分块矩阵的求逆公式如下(6)：设R矩阵是一个非奇异矩阵， R-〔8H, 则 E-1+E-1FDGE-1-E-1FD-1 R-1= -D-1GE-1 D-1 其中D=H-GE~1F。把(9)式代入此公式，并设 R。1=(XXn)1,有 R:1+-1-RiiXi1ZaZiX.Ri1--1-R-1nXiZo R1a+1=(XTn+1Xm+1)1= Ca an -1ZiXRal 1 (11) 其中aa=ZZn-ZX.R1XZn。 (12) 这就建立了AR(n)和AR(n+1)模型之间的递推关系。由求AR(n)模型时已知 的Ra1=(XXm)~1递推求出(X"n+1Xm+1)-1,进而求得AR（n+1)模型，从而减少了 计算量。目前一般的AR模型建模策略是：先逐个求出AR(1)、AR(2)、…AR (如ax)这样若千个模型，然后根据一定的检验准则，并考虑其它因素，如AR谱的质量，模型 的精度等，最后选定其中某一阶模型AR(n)(1≤n≤nmax),作为最终模型。本文提出 的方法适用于这二过程。取预留数据个数m≥nmax即可。 除样本很短，没有预留数据的余地外，该方法是可行的。虽然当模型阶数变化时，数 据的取用长度有变化，但变化很小。同时就每一阶模型而言，无论等于几，都是独立地由 最小二乘估计得出的，在残差平方和最小的意义下都是精确的，没有损失估计精度。 2减少计算量的其它措施 我们研究了各运算过程量的多种形式、规律和特点，以及以往常用的最小二乘法计算机 程序，认为有下述一些关系和方法可用来进一步减少计算工作量。 (1)XTX矩阵是对称矩阵，由线性代数的基本定律可知R-1=(XTX)·1也是对称矩 阵，因此只需计算它的上三角区元素就可以了。 (2)由于采用递补数据的办法，有一系列递推关系可以利用，这是本文方法突出的优 点。 设 r(0,0) (0;1)…r(0,n-1) XIX= r(1,0) r(1,1)…r(1,n-1) (13) r(n-1,0)r(n-1,1)…r(n-1,n-1) N1 其中r（i,j)=∑X-iX-j。 (14) t=0 51
矩 阵的逆 万 。 一 ‘ 构 成递推关 系 。 式 则直接给 出了 。 与 二 的关系 。 厂仗 、尸 分 块粗阵的求 逆 公式如下 一 〔 〕 设 矩 阵 是一个非 奇 异矩 阵 ， 则 一 〔 一 一 一 一 一 尸 一 ’ 〕 其中 一 一 ‘ 。 把 式代 入此公式 ， 并设 石 二 百 。 一 ‘ ， 有 一 。 一 一 ， ， 百 。 孟 一 一 一 。 万 。 一鱼一 王 孟‘ 其中 。 百 。 一 蕊 ‘ 石‘ 百 这 就建立 了 和 十 模型 之 间的递 推关 系 。 由求 模 型 时 已知 的 孟‘ 王 。 一 ‘ 递 推求 出 、 一 ‘ ， 进而 求 得 模型 ， 从而减少 了 计算量 。 目前一般的 模型 建模 策 略 是 先 逐 个 求 出 、 、 … 二 这 样若 干个模型 ， 然后根据一 定 的检验准 则 ， 并考虑 其它 因素 ， 如 谱 的质量 ， 模型 的精度等 ， 最后选定 其 中某一 阶模型 《 《 ， 作为最终模型 。 本文提 出 的方法适 用于这一 过 程 。 取预 留数据个数 二 即可 。 除样本很短 ， 没 有预 留数 据 的 余地 外 ， 该 方法 是可行 的 。 虽然 当模型 阶数 变化时 ， 数 据的取用长度有 变化 ， 但 变化 很小 。 同时就每一 阶模型而言 ， 无论 等于几 ， 都 是独 立地 由 最小 二乘估计得 出的 ， 在 残差 平方和 最小 的意义 下都是精确 的 ， 没 有损失估计精度 。 减少计算量 的其它措施 我 们研究 了各运算过程量 的 多种 形 式 、 规律和特 点 ， 以 及以 往常 用 的最小二乘法计算机 程 序 ， 认为有下述 一些关 系和方法 可 用来进一步减 少计算工 作量 。 ’ 矩 阵是对称矩 阵 ， 由线 性代数 的基本定律可知 一 ‘ ’ 一 ‘ 也 是对称矩 阵 ， 因此 只需 计算它 的上三角区元 素就可 以 了 。 由于采 用递 补数据的办法 ， 有一 系列递推关 系可 以利 用 ， 这 是 本文方法突 出的优 点 。 万 。 ， 万 … ， ， … 一 ， 一 ， … 一 ， 一 其中 ， 一 艺 一 一
又 r(0,0)…r(0,n-1)1r(0,n) X+1X+1= r(n-1,0)…r(n-1,n-1)1r(n-1,n) 1---- 〔经. r n,0)...r n,n-1)I r(n,n) N一1 于是ZZ。=r（n,n)=∑X2t-,考虑到由XTa-1Xm-1递推到XX时有Zn-1Zm-1 t0 N-1 =r（n-1,n-1)=∑X2t-n+1。有递推关系 t=0 N-1 N-1 Z日Z。=zX2-n=ΣXt-n+1-X2N-n+X2_n=ZTn-1Zm-1-XN-n+X:。 t=0 t=0 (16) 由(12)武有XhY=〔Y门. N-1 每次升阶时只需计算ZY=∑Xt+1Xt- t=0 =r（-1,n)。这个值必须由样本直接计算。 类似(16)，对于XZ有递推关系 N-1 r(0,n) X,Xt-n t=0 -Xx+XoX t=0 N-1 N-1 XIZ=r(1,n) Xt-Xt-n = XtXt-n+1-XN-1XN-+X_iX. t=0 t=0 t 年 N-1 N-1Xt-n+2X-n+1-XN-n+1XN-n+ r(n-1,n) X-+X- X-n+1X-a r(-1,n-1) XN Xo r(0,n-1) XN- XN-1 +X-n X-1 (r(n-2,n-1) XN-+1 X-n+1 ZTn-1Y XN Xo XN-1 X一1 XTn-1Zn-1 XN-n +X-n (17) XN-+1 X-n+1) (3)虽然每次升阶时只需要由样本数据直接计算ZY这一个数，但计算量是比较多的， 尤其当数据量N较大时更甚。这里引入了Winograd向量内积快速算法()，用它来计算多个 向量的内积时，大约可以减少一半的乘法运算量。由下边的推导可以看出，引入这一算法的 收益是很大的。Winograd:算法的原理如下： 52
又 ， … ， 一 ， 、 十 、 十 一 ， 一 一 ， 一 一 ， 荟 。 艺 。 盖 。 万 。 〕 ， … ， 一 ， 于 是 万 。 ， 二 艺 “ 卜 。 。 考虑到由 一 一 递推到 荟 。 时有 ， 一 一 二 一 ， 一 二 艺 “ 卜 。 。 有递推关系 万 。 。 一。 一 一。 。 。 一 一 一 一 兰 。 由 式有 十 每 次升阶时只需计算 王 二 一 。 一 ， 。 这个值必须 由样本直接计算 。 类似 ， 对于 王 有递推关 系 一… ， 一 一 一 ， 一 。 。 ， 一 艺 一 一 一 一 一 一 。 一 。 十 一。 ‘ 、苦分 一 、 一 一 。 一 。 ， ， 一 一 一。 ‘ 一一 、 ‘ 产卫‘ 一 ， 一 ， 一 ， 一 一 ， 一 一… 一 一 一 一 一 。 一 。 丁 。 一 。 一 一 一 一 一 一 一 一 虽 然每 次升 阶 时 只需 要 由样本数据直接讨算 工 这 一个数 ， 但计算量 是 比较 多的 ， 尤 其 当数 据量 较大时 更甚 。 这 里 引人 了 向量 内积快速算法〔 〕 ， 用它 来计算多个 向量的 内积时 ， 大 约可 以减少一半的乘 法运算量 。 由下边 的推导可 以看 出 ， 引入这一算法 的 收 益 是很大的 。 算法的原理如下
设有两个向量 a 1 A B、 内积 ATB=ab1+a2b2+…+anbw 先计算A向量和B向量的奇偶积门A和I, CN/2) cN12), a21-1a21y j=1 n，b-b2。 j=1 其中记号(N/2)={N/2 N为偶数 (N-1)/2 N为奇数。 〔N/2) (a2j-1+b2i)（azi+b2j-1)-na-n。 有ATB= j=1 N为偶数 (N/2) (j=I a2j-1+b2j)(a2j+b2j-1)-nA-n.+aNbN N为奇数 CN/2) 令0= ∑（a:-1+b2i)(a2j+b2j-1) j=1 (18) ·.AB={0_二i+aba N为偶数 N为奇数 (19) 在由n=1,2,,nmax逐阶递推建模的全过程中，共要由样本直接计算X1X1、 XTY和ZY(n=1,…,nma1这nmax+1个内积。其中 Xo X:= X-1 ,Y= :,Z:= …yZ。 XN-1 XN XN-2 XN--1 当N为偶数时，设Y、X:、乙n的奇偶积分别为n￥、门x、n,则： ny=X:X2+XaX4+...+Xx-:Xx> (20) 门x=XoX1+X2X3+…+Xw-2XN-1, (21) 存在如下递推关系 n1=门r-XN-1Xw+X.1Xo n2=nx-XN-2XN-1+X_2X_1 (22) nn=nn-2-XN-nXN-+1+X-nX-n+1 然后利用(18)、（19)两式就可求出所需的nmx+1个向量积，这里不再赘述。这样 共用到Nn/2+8N/2+2n-2次乘法和8Nn/2+5N/2+3n-3次加法运算。若直 接计算要用Nn+N次乘法和Nn+N-n-1次加法。乘法运算减少了将近-一半，加法计算虽 有增加，但对于计算机来说，做加法比做乘法运算快得多。这正是Winograd算法的特点。 (4)以上给出了绝大部分计算公式。在建模工作中往往还要估计残差平方和E。。有 53
设有两个向量 二 … 三 ， 三 贬 一 二 气 内积 … 、 先计算 向量和 向量的奇偶积” 人 和” ， 〔 幻 刀 艺 卜 一， 其 中记号 〔 〕 气 了 一 〔 〕 ” 。 艺 一 ， 为偶数 为奇数 。 有 〕 艺 一 通一 一 刀 一 叮 〔 〕 艺 一 十 一 一 ” 一 ” 。 为偶数 为奇数 〔 〕 令 。 艺 一 卜， 。 自 一 一 刀 一 月 一 刀 。 为偶数 为奇数 在 由 ， ， … ， 二 逐阶递推建模的全过程 中 ， 共要 由样本直 接 计 算 ， 、 和 五 ， … ， 。 ‘ 这 。 。 十 个 内积 其中 。 ， ‘ “ 三 卜 “ 三 卜 ‘ 气 一 夕 少 ， ， ” ， · 一 。 一 一 当 为偶数时 ， 设 、 、 的奇偶积分别为 丫 、 。 、 。 ， 则 丫 … ， ， ” … 、 卜 、 ， 存在如下递推关系 厂 刀 ， 勺， 一 一 月 月 一 一 一 刀 。 刀。 一 一 一 一 一 一 然后利 用 、 两式就可求 出所需的 个 向量积 ， 这里 不再赘述 。 这样 共用到 十 一 次乘法 和 一 次 加法运算 。 若 直 接计算要用 十 次乘法和 十 一 一 次加 法 。 乘法运算减 少 了将近 一半 ， 加法计算虽 有 增加 ， 但对于计算机来说 ， 做加法比做乘法运算快得 多 。 这正 是 算法 的特 点 。 以 上给 出了绝大部分计算公式 。 在建模工 作 中往往还要估计残差平方和 。 。 有
仓=芝a=Y-可：XY (23) t=1 此式中YtY与n无关，只需计算一次。本。，X:Y都在建模中得出，因此很容易求出E值。 这种递推求逆建立各阶AR模型的方法，包括求出各残差平方和估计值仓，1，R（n=1, 2…,n)共用到Nn/2+2N+5n3/6+5n2/2+11n/3-5次乘法和3Nn/2+4N +5n3/6+n2+19n/6-3次加法运算。 3 第二种方法 预留少量数据递补进入计算，使X'a+1Xn+1阵的左上角n×n阶阵等于X日Xn,建立了 AR(n+1)和AR(n)之间的关系。由此还可以导出更为简捷的建模方法。。=（X: Xn)ˉ1XY式的原形为XXn中。=XY是一个线性方程组。Cout分解法是解线性方程组的 一种常用方法，它的原理是(8)： a 11..a in 设AX=b,其中A= X- aa1…aan X b。 把A分解为L、U两个三角阵，即 1**41n A=LU,L U= la1…lanf 0 .LUX=b。 令UX=y,y=〔y,…ya, (24) .Ly=b (25) 由.(25)式进代求出中间量y,再由(24)式迭代求出未知数X。在把A矩阵分解为L、U两 个三角阵时有公式 j-1 1)=a11-Σ.11xuxj i≥j k=1 i-1 u11=a11-,Σlii“Ki K<J 上=1 (u:i=u1i=1 这三个公式结合(9)式不难证明，在XX。=An=L.Ua,X7n+1Xn+1=An+1= Ln+1Un+1时，不仅An+1的左上角n×n阶方阵就是Aa,L+1和Un+1的左上角n×n阶阵也 分别是L。和Um。即 54
洲、 、 付 尹、 、 。 ， ， 、 艺 蓄二 一 巾 工 王 此式中 · 与 无关 ， 只需计算一次 。 多 。 、 都在建模中得出 ， 因此很容易求 出成值 。 这种递 推求 逆建立 各阶 模型的方法 ， 包括求 出各残差平方和估计 值 。 ， 、 ， 、 ， … ， 共用到 “ 一 次 乘 法 和 。 “ 一 次加法运算 。 第二种方法 预 留少量数据递 补进 人计算 ， 使 几十， 十 阵的左上 角 阶 阵等于 二 。 ， 建立 了 和 之 间的关系 。 由此还 可 以导 出更为简捷的建模方 法 。 巾 。 盖 一 ’ 百 式的原形 为 二 公 。 二 百 是一个线性方程组 。 分解法 是解线性方程组 的 一 种常用方法 ， 它 的原理 是 ⑧ 玩 · 饥 一 … 」 设 一 ， 其 中 二 三 二 ‘ 。 把 分解为 、 两个三 角阵 ， 即 ︸、 。 ，… 。 。 ，老 一一 一 。 令 ， 〔 ，， … 。 〕 ， 二 式迭代求 出未知 数 。 在把 矩 阵分解为 、 两 ， 。 ，，一 】 一 ， ‘ 遥 “ 》 由 式迭代求 出 中间量， 酗向问归、。 个 三角阵时有 公 式 盆 五 这 三个 公式结合 式不难 证 明 ， 在 万 。 二 。 十 时 ， 不仅 。 十 的左 上角 火 阶方阵就 是 。 ， 分别 是 。 和 。 即 。 。 ， 。 士 。 和 。 的左上角 阶阵 也
1u1n-+1 e L UI: Ln+1= I un,n+i 0…0 1 这样在求AR(n+1)的系数中a+1时，可以利用已知的Ln、Un,只需增加计算(ln+11, …,ln+1n十1)和(u1n+1,…,unn+1)这些数这使计算工{，i、.求逆的办法进一步减小。 在其它计算仍采用前述方法的情况下，完成AR(1),,AR(n)的建模，包括求出各 阶的En共需要Nn/2+2N+n3/2+5n/2+4n-4次乘法和3Nn/2+4N+n3/2 +2n2+7n/2一3次加法运算。 这个方法还具有其它优点：它的计算机程序简单，由于它的算法简捷，因而引入的计算 误差小；它的L、Un两矩阵实际上包含了各阶的L:、U:(i≤n)矩阵，由这两个矩阵可以 直接求出我们指定的任何一阶AR（i)(i≤n)的模型系数，这会带来一定的方便。 另外，还可以采用平方根法、Cholesky分解法等其它算法来求解。它们此Crout分解法 计算量稍大，但由于XTX阵是对称的，相应地可减少部分所占用的计算机内存。在需要时可 以采用。·。 4 实用效果 本文给出的两种方法的计算量与Burg算法、Marple算法计算量的比较见表。表中 Burg算法按3Nn-n2-N+3n次乘法和3Nn-n-2N-n次加法运算，Marple算法按 Nn+9n2+2N+25n-3次乘法和Nn+8n2+N+7n-8次加运算计算()。普通最小 二乘法完成前述计算要用到比例于N3+n4次的运算。在利用了各计算量之间的关系加以改 进后（如美国威斯康星大学1982年左右的程序），仍需比例于Nn+n4次的运算。与之相 比，本文方法减少了约为倍（几倍到几十倍）的计算工作量。 计算量的比较 Table Comparison of computational operations Methods Burg/a Marple/s First Second n multiplication addition 10 100 1:0.83:0,64:0,53 1:0,73:1.06:0.98 20 1:1.13:1.61:1,13 1:1.01:1.96:1.54 10 250 1:0.58:0.40:0.35 1:0,52:0,83:0,80 20 1:0.67:0.740.56 1:0,61:1.11:0,95 500 10 1:0,50:0,32:0.30 1:0.46:0,75:0.74 30 1:0.57:0,76:0.56 1:0.53:1,11:0.93 1000 10 1:0.45:0,28:0.27 1:0.43:0,72:0,71 30 1:0.46:0,47:0.37 1:0.44:0.83:0,74 55
、 … 了 崖挤 一卜 ， 十 “一 …一 、 一 ， 十 下卜 产 ‘ “ ， 一 几一 丁一 仁 一 丁 一 ‘ ” ， ” ‘ ’ ，” 一 ‘ ‘ 二 这 样在求 十 的 系数中 十 时 ， 可以利 用 已知 的不 。 、 ， 只需增加 计算 ， ， … ， ，叶 ‘ 和 ， ， … ， “ 。 这 些数 。 这使计算工 二 、 求逆 的办法进 一步减小 。 在 其它计算仍采 用前述 方法的情况下 ， 完成 ， 。 二 ， 入坟 的建模 ， 包 括求 出各 阶的宜 。 共需要 。 十 十 。 。 、 一 次乘法和 。 十 。 “ 十 一 次加法运算 。 这个方法还具有 其它优点 它 的计算机程序简单， 由于它的算法简捷 ， 因而引人 的计算 误 差小 它 的 、 两矩 阵实际上包含了各阶的 、 矩 阵 ， 由这 两个矩 阵可 以 直接求 出我们指定的 任何，阶 ‘ 《 的模型 系数 ， 这 会带 来一定 的方便 。 另外 ， 还可 以采 用平方根法 、 分解法等 其它算法 来求解 。 它 们 比 分解 法 计算量稍大 ， 但由于 阵是对称的 ， 相应地可减少部分 所 占用的计算机 内存 。 在需要时可 以采 用 。 实 用 效 果 本文给 出的两种 方法 的计算量与 算法 、 算法计算量 的 比 较 见 表 。 表 中 算法 按 一 “ 一 次乘法 和 一 “ 一 一 次加法运算 ， 算 法 按 一 次乘法和 “ ，卜 一 次加运 算计算〔 〕 。 普 通 最 小 二乘法 完成前述计算要用到 比例 于 “ 次的运算 。 在 利 用 了各计算量 之 间的关 系加 以改 进后 如美国威斯康星大学 年左 右的程 序 ， 仍需比例 于 “ 十 毛次 的运算 。 与 之相 比 ， 本文方法减 少 了约为 倍 几 倍到几 十倍 的计算工 作量 。 计 算 量 的 比 较 ， 三 ， 一 。 名 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。
这两种方法均已编成计算程序，使用效果良好。建模的我差平方和E。与Marple算法 的结果很接近，.有些情况下优于后者。所得到的AR谱与Marple算法的AR谱比较，未见明 显差别。本文提出的两种方法的结果十分接近，只有因引入的计算误差不同造成的微小差 别。另外，由表可以看出，这两种方法在样本数据量大的情况下优点更突出。 参考文献 1 Pandit,S.M.;Wu,S.M.Time Series and System Analysis'with Applications )John Wiley and Sons,New York,1983 〔2〕安鸿志等：《时间序列的分析与应用》，科学出版社，1983 2 3 )Kay S.M.Marple S.L.Spectrum Analysis-A Modern Perspective, Proc.IEEE,Vol.69(1981),1380 [4 Morf M.Vieira A.;Lee D.T.;Ladder Forms for Idenfication and Speech Pocessing,.in Proc.1977,Conf.Decision and Control,P.1074- 1078,Dec.1977 [5 MarpleL.:A New Autoregressive Spectrum Analysis Algorithm,IEEE Irans,Acon,Speech,Sig.Proc Vol,Assp-28,No.4 P.441-453,Aug.1980 〔6)夏天长：系统辨识，清华大学出版社，1983,34 〔7〕游兆永：线性代数与多项式的快速算法，上海科学技术出版社，1980,1一5 〔8〕华中工学院数学教研室：算法语言。计算方法，人民教育出版社，1978,231 56
这两 种 方法 均 已编 成计算程 序 ， 使 用效果 良好 。 建模的残差平方 和雹 与 。 算 法 的结果很接近 ， 有些 情况下优 于后者 。 所得 到 的 谱与 算法 的 谱比较 ， 未见明 显差别 。 本文提 出的两 种方法的结果十分接近 ， 只有因引人 的计算误 差不 同造成 的微 小差 别 。 另外 ， 由表可 以看 出 ， 这币种方法在样本数据量大 的情况下优 点更突 出 。 〔 〕 ， ， 参 考 文 献 《 ‘ 》 ， ， ， 〔 〕 安鸿志 等 《 时 间序列 的分析与应用 》 ， 科学 出版社 ， 〔 〕 挂 一 “ ， 飞 ， ， 〔 〕 ， ， ， 〔 〕 ， ， ， ， ， ， 一 一 ， 一 ， 〔 〕 夏天长 系统辨识 ， 清华大 学 出版社 ， ， 〔 〕 游 兆 永 线 性代数与 多项式的快速算法 ， 上海科学 技术 出版社 ， ， 一 〔 〕 华中工 学院数学 教研室 算法语言 计算方法 ， 人 民教育 出版社 ，
