A Levinson solver for symmetric square Toeplitz systems of Complex
type.
This class provides members for inverting a symmetric square Toeplitz matrix
(see ComplexDoubleSymmetricLevinson.GetInverse member), calculating the determinant of the matrix
(see ComplexDoubleSymmetricLevinson.GetDeterminant member) and solving linear systems associated
with the matrix (see Solve
members).
The class implements a UDL decomposition of the inverse of the Toeplitz matrix. The decomposition is based upon Levinson's algorithm. As a consequence, all operations require approximately N squared FLOPS, where N is the matrix order. This is significantly faster than Cholesky's factorization for symmetric matrices, which requires N cubed FLOPS.
A requirement of Levinson's algorithm is that all the leading sub-matrices and the principal matrix must be non-singular. During the decomposition, sub-matrices and the principal matrix are checked to ensure that they are non-singular. When a singular matrix is found, the decomposition is halted and an internal flag is set. The ComplexDoubleSymmetricLevinson.IsSingular property may be used to access the flag, to determine if any singular matrices were detected.
A outline of this approach to the UDL decomposition of inverse Toeplitz matrices is found in the following reference:
Sun-Yuan Kung and Yu Hen Hu, A Highly Concurrent Algorithm and Pipelined Architecture for Solving Toeplitz Systems, IEEE Transactions on Acoustics, Speech and Signal Processing, Volume ASSP-31, Number 1, Febuary 1983, pages 66 - 75.
Copyright (c) 2003-2004, dnAnalytics Project. All rights reserved. See http://www.dnAnalytics.net for details.
Adopted to Altaxo (c) 2005 Dr. Dirk Lellinger.
The following simple example illustrates the use of the class:
using System; using dnA.Exceptions; using dnA.Math; using System.IO; namespace Example_5 { class Application { // application entry point [STAThread] public static void Main(string[] args) { ComplexDoubleVector cdv = new ComplexDoubleVector(4); cdv[0] = new Complex(4.0, 0.0); cdv[1] = new Complex(12.0, -4.0/3.0); cdv[2] = new Complex(80.0/3.0, 64.0/3.0); cdv[3] = new Complex(48.0, -16.0/3.0); // create Levinson solver ComplexDoubleSymmetricLevinson cdsl = new ComplexDoubleSymmetricLevinson(cdv); // display the Toeplitz matrix ComplexDoubleMatrix T = cdsl.GetMatrix(); Console.WriteLine("Matrix:: {0} ", T.ToString("E3")); Console.WriteLine(); // check if matrix is singular Console.WriteLine("Singular: {0}", cdsl.IsSingular); // get the determinant Console.WriteLine("Determinant: {0}", cdsl.GetDeterminant().ToString("E3")); Console.WriteLine(); // get the inverse ComplexDoubleMatrix Inv = cdsl.GetInverse(); Console.WriteLine("Inverse:: {0} ", Inv.ToString("E3")); Console.WriteLine(); // solve a linear system ComplexDoubleVector Y = new ComplexDoubleVector(4); Y[0] = new Complex(1036.0/3.0, -212.0); Y[1] = new Complex(728.0/3.0, -200.0/3.0); Y[2] = new Complex(388.0/3.0, -148.0/3.0); Y[3] = new Complex(304.0/3.0, -40.0/3.0); ComplexDoubleVector X = cdsl.Solve(Y); Console.WriteLine("X:: {0} ", X.ToString("E3")); Console.WriteLine(); } } }
The application generates the following results:
Matrix:: rows: 4, cols: 4 4.000E+000 + 0.000E+000i, 1.200E+001 -1.333E+000i, 2.667E+001 + 2.133E+001i, 4.800E+001 -5.333E+000i 1.200E+001 -1.333E+000i, 4.000E+000 + 0.000E+000i, 1.200E+001 -1.333E+000i, 2.667E+001 + 2.133E+001i 2.667E+001 + 2.133E+001i, 1.200E+001 -1.333E+000i, 4.000E+000 + 0.000E+000i, 1.200E+001 -1.333E+000i 4.800E+001 -5.333E+000i, 2.667E+001 + 2.133E+001i, 1.200E+001 -1.333E+000i, 4.000E+000 + 0.000E+000i Singular: False Determinant: -1.478E+006 - 8.879E+005i Inverse:: rows: 4, cols: 4 -3.427E-003 + 2.015E-003i, 3.534E-003 + 2.723E-003i, 1.145E-002 - 2.307E-002i, 2.735E-003 + 7.280E-003i 3.534E-003 + 2.723E-003i, -1.702E-002 + 2.724E-005i, 1.534E-002 + 3.027E-002i, 1.145E-002 - 2.307E-002i 1.145E-002 -2.307E-002i, 1.534E-002 + 3.027E-002i, -1.702E-002 + 2.724E-005i, 3.534E-003 + 2.723E-003i 2.735E-003 + 7.280E-003i, 1.145E-002 - 2.307E-002i, 3.534E-003 + 2.723E-003i, -3.427E-003 + 2.015E-003i X:: Length: 4 1.000E+000 - 1.000E+000i, 2.000E+000 + 2.000E+000i, 3.000E+000 - 3.000E+000i, 4.000E+000 - 4.000E+000i
Constructor | Description | ||||
Full Usage:
ComplexDoubleSymmetricLevinson(T)
Parameters:
IROComplexDoubleVector
-
The left-most column of the Toeplitz matrix.
|
Constructor with
|
Instance member | Description | ||||||
|
Get the diagonal matrix of the UDL factorisation. It is recommended that the ComplexDoubleSymmetricLevinson.IsSingular property be checked to see if the decomposition was completed, before attempting to obtain the diagonal matrix.
|
||||||
|
Get the determinant of the Toeplitz matrix. It is recommended that the ComplexDoubleSymmetricLevinson.IsSingular property be checked to see if the decomposition was completed, before attempting to obtain the determinant.
|
||||||
|
Get the inverse of the Toeplitz matrix. The class implicitly decomposes the inverse Toeplitz matrix into a UDL factorisation using the Levinson algorithm, before using Trench's algorithm to complete the calculation of the inverse. Trench's algorithm requires approximately N squared FLOPS, compared to N cubed FLOPS if we simply multiplied the UDL factors (N is the matrix order).
|
||||||
|
Get a copy of the Toeplitz matrix.
|
||||||
|
Get a vector that represents the left-most column of the Toplitz matrix.
|
||||||
Full Usage:
this.IsSingular
Returns: bool
|
Check if the Toeplitz matrix or any leading sub-matrices are singular. If the Toeplitz matrix or any leading sub-matrices are singular, it is not possible to complete the UDL decomposition using the Levinson algorithm.
|
||||||
|
Get the lower triangle matrix of the UDL factorisation. It is recommended that the ComplexDoubleSymmetricLevinson.IsSingular property be checked to see if the decomposition was completed, before attempting to obtain the lower triangle matrix.
|
||||||
Full Usage:
this.Order
Returns: int
|
Get the order of the Toeplitz matrix.
|
||||||
Full Usage:
this.Solve
Parameters:
IROComplexDoubleVector
-
The right-hand side of the system.
Returns: ComplexDoubleVector
The solution vector.
|
Solve a symmetric square Toeplitz system with a right-side vector. This member solves the linear system TX = Y, where T is the symmetric square Toeplitz matrix, X is the unknown solution vector and Y is a known vector. The class implicitly decomposes the inverse Toeplitz matrix into a UDL factorisation using the Levinson algorithm, and then calculates the solution vector.
|
||||||
Full Usage:
this.Solve
Parameters:
IROComplexDoubleMatrix
-
The right-hand side of the system.
Returns: ComplexDoubleMatrix
The solution matrix.
|
Solve a symmetric square Toeplitz system with a right-side matrix. This member solves the linear system TX = Y, where T is a symmetric square Toeplitz matrix, X is the unknown solution matrix and Y is a known matrix. The class implicitly decomposes the inverse Toeplitz matrix into a UDL factorisation using the Levinson algorithm, and then calculates the solution matrix.
|
||||||
|
Get the upper triangle matrix of the UDL factorisation. It is recommended that the ComplexDoubleSymmetricLevinson.IsSingular property be checked to see if the decomposition was completed, before attempting to obtain the upper triangle matrix.
|
Static member | Description | ||||||
Full Usage:
ComplexDoubleSymmetricLevinson.Inverse(T)
Parameters:
IROComplexDoubleVector
-
The left-most column of the symmetric Toeplitz matrix.
Returns: ComplexDoubleMatrix
The inverse matrix.
|
Invert a symmetric square Toeplitz matrix. This static member combines the UDL decomposition and Trench's algorithm into a single algorithm. When compared to the non-static member it requires minimal data storage and suffers from no speed penalty. Trench's algorithm requires N squared FLOPS, compared to N cubed FLOPS if we simply solved a linear Toeplitz system with a right-side identity matrix (N is the matrix order).
|
||||||
Full Usage:
ComplexDoubleSymmetricLevinson.Solve(T, Y)
Parameters:
IROComplexDoubleVector
-
The left-most column of the Toeplitz matrix.
Y : IROComplexDoubleVector
-
The right-side vector of the system.
Returns: ComplexDoubleVector
The solution vector.
|
Solve a symmetric square Toeplitz system with a right-side vector. This method solves the linear system AX = Y. Where T is a symmetric square Toeplitz matrix, X is an unknown vector and Y is a known vector. This static member combines the UDL decomposition and the calculation of the solution into a single algorithm. When compared to the non-static member it requires minimal data storage and suffers from no speed penalty.
|
||||||
Full Usage:
ComplexDoubleSymmetricLevinson.Solve(T, Y)
Parameters:
IROComplexDoubleVector
-
The left-most column of the Toeplitz matrix.
Y : IROComplexDoubleMatrix
-
The right-side matrix of the system.
Returns: ComplexDoubleMatrix
The solution matrix.
|
Solve a symmetric square Toeplitz system with a right-side matrix. This method solves the linear system AX = Y. Where T is a symmetric square Toeplitz matrix, X is an unknown matrix and Y is a known matrix. This static member combines the UDL decomposition and the calculation of the solution into a single algorithm. When compared to the non-static member it requires minimal data storage and suffers from no speed penalty.
|
||||||
Full Usage:
ComplexDoubleSymmetricLevinson.YuleWalker(R)
Parameters:
IROComplexDoubleVector
-
The left-most column of the Toeplitz matrix.
Returns: ComplexDoubleVector
The solution vector.
|
Solve the Yule-Walker equations for a symmetric square Toeplitz system This member is used to solve the Yule-Walker system AX = -a, where A is a symmetric square Toeplitz matrix, constructed from the elements R[0], ..., R[N-2] and the vector a is constructed from the elements R[1], ..., R[N-1]. Durbin's algorithm is used to solve the linear system. It requires approximately the N squared FLOPS to calculate the solution (N is the matrix order).
|