FloatSymmetricLevinson Type
A Levinson solver for symmetric square Toeplitz systems of float type.
This class provides members for inverting a symmetric square Toeplitz matrix
(see FloatSymmetricLevinson.GetInverse member), calculating the determinant of the matrix
(see FloatSymmetricLevinson.GetDeterminant member) and solving linear systems associated
with the matrix (see Solve members).
The class implements an 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 FloatSymmetricLevinson.IsSingular property may be used to access the flag, to determine if any singular matrices were detected.
It has been shown that Levinson's algorithm is weakly numerically stable for positive-definite Toeplitz matrices. It is usual to restrict the use of the algorithm to such matrix types. The FloatSymmetricLevinson.IsPositiveDefinite property may be checked to verify that the Toeplitz matrix is positive-definite.
If one of the leading sub-matrices or the principal matrix is near-singular, then the accuracy of the decomposition will be degraded. An estimate of the resulting error is provided and may be accessed with the FloatSymmetricLevinson.Error property. This estimate is only valid for positive-definite matrices.
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.
Example
The following simple example illustrates the use of the class:
using System;
using dnA.Exceptions;
using dnA.Math;
using System.IO;
namespace Example_3
{
class Application
{
// format string for matrix/vector elements
private const string frmStr = " 0.000E+000;-0.000E+000";
// application entry point
[STAThread]
public static void Main(string[] args)
{
// create Levinson solver
FloatSymmetricLevinson fsl = new FloatSymmetricLevinson(1.0f, 0.5f, 0.2f);
// display the Toeplitz matrix
FloatMatrix T = fsl.GetMatrix();
Console.WriteLine("Matrix:: {0} ", T.ToString(frmStr));
Console.WriteLine();
// check if matrix is singular
Console.WriteLine("Singular: {0}", fsl.IsSingular);
// check if matrix is positive definite
Console.WriteLine("Positive Definite: {0}", fsl.IsPositiveDefinite);
// get error for inverse
Console.WriteLine("Cybenko Error Estimate: {0}", fsl.Error.ToString("E3"));
// get the determinant
Console.WriteLine("Determinant: {0}", fsl.GetDeterminant().ToString("E3"));
Console.WriteLine();
// get the inverse
FloatMatrix Inv = fsl.GetInverse();
Console.WriteLine("Inverse:: {0} ", Inv.ToString(frmStr));
Console.WriteLine();
// solve a linear system
FloatVector X = fsl.Solve(4.0f, -1.0f, 3.0f);
Console.WriteLine("X:: {0} ", X.ToString(frmStr));
Console.WriteLine();
}
}
}The application generates the following results:
Matrix:: rows: 3, cols: 3
1.000E+000, 5.000E-001, 2.000E-001
5.000E-001, 1.000E+000, 5.000E-001
2.000E-001, 5.000E-001, 1.000E+000
Singular: False
Positive Definite: True
Cybenko Error Estimate: 4.087E-007
Determinant: 5.600E-001
Inverse:: rows: 3, cols: 3
1.339E+000, -7.143E-001, 8.929E-002
-7.143E-001, 1.714E+000, -7.143E-001
8.929E-002, -7.143E-001, 1.339E+000
X:: Length: 3
6.339E+000, -6.714E+000, 5.089E+000
Constructors
| Constructor | Description | ||||
Full Usage:
FloatSymmetricLevinson(T)
Parameters:
IReadOnlyList<float32>
-
The left-most column of the Toeplitz matrix.
|
Constructor with
|
||||
Full Usage:
FloatSymmetricLevinson(T)
Parameters:
float32[]
-
The left-most column of the Toeplitz matrix.
|
Constructor with
|
Instance members
| Instance member | Description | ||||||
|
Get the diagonal matrix of the UDL factorisation. It is recommended that the FloatSymmetricLevinson.IsSingular property be checked to see if the decomposition was completed, before attempting to obtain the diagonal matrix.
|
||||||
Full Usage:
this.Error
Returns: float32
|
Get an error estimate for the inverse matrix. This estimate is approximate and is only valid for positive definite matrices. It is useful for checking the accurary of the UDL decomposition. If the Toeplitz matrix or one of the leading sub-matrices is near-singular, the accuracy of the decomposition will be degraded. This error estimate can be used to identify such situations. The estimate is based on the Cybenko error-bound for the inverse and is a relative error. When calculating this error it is assumed that the Single type utilises a 23-bit mantissa. The error estimate will range in value between 1.19e-7 and Single.PositiveInfinity.
ExampleThe following example calculates the inverse of a symmetric positive-definite Toeplitz matrix. An exact inverse was obtained previously using a computer algebra system. The exact inverse is compared to the inverse calculated by this class. The actual error is identified and compared to the Cybenko error-bound.
using System;
using dnA.Exceptions;
using dnA.Math;
using System.IO;
namespace Example_4
{
class Application
{
// The main entry point for the application.
[STAThread]
public static void Main(string[] args)
{
// exact inverse from yacas
FloatMatrix exact = new FloatMatrix(5);
exact[4, 4] = exact[0, 0] = +1.3577842E+000f;
exact[3, 4] = exact[4, 3] = exact[1, 0] = exact[0, 1] = -5.9039646E-001f;
exact[2, 4] = exact[4, 2] = exact[2, 0] = exact[0, 2] = -1.0830325E-001f;
exact[1, 4] = exact[4, 1] = exact[3, 0] = exact[0, 3] = -5.9423022E-002f;
exact[4, 0] = exact[0, 4] = -5.8145106E-002f;
exact[3, 3] = exact[1, 1] = +1.6120124E+000f;
exact[3, 2] = exact[2, 3] = exact[2, 1] = exact[1, 2] = -5.4584837E-001f;
exact[3, 1] = exact[1, 3] = -8.7102652E-002f;
exact[2, 2] = +1.6180506E+000f;
// create levinson solver
FloatSymmetricLevinson fsl = new FloatSymmetricLevinson(1.0f, 1.0f/2.0f, 1.0f/3.0f, 1.0f/4.0f, 1.0f/5.0f);
// identify relative error
FloatMatrix Error = fsl.GetInverse() - exact;
float e = Error.GetL1Norm() / exact.GetL1Norm();
// display results
Console.WriteLine("Observed Error: {0}", e.ToString("E3"));
Console.WriteLine("Cybenko Error Estimate: {0}", fsl.Error.ToString("E3"));
Console.WriteLine();
}
}
}The application generates the following results.
Observed Error: 6.110E-008
Cybenko Error Estimate: 5.520E-007
|
||||||
Full Usage:
this.GetDeterminant
Returns: float32
The determinant.
|
Get the determinant of the Toeplitz matrix. It is recommended that the FloatSymmetricLevinson.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.IsPositiveDefinite
Returns: bool
|
Check if the Toeplitz matrix is positive definite. It has only been shown that the Levinson algorithm is weakly numerically stable for symmetric positive-definite Toeplitz matrices. Based on empirical results, it appears that the Levinson algorithm gives reasonable accuracy for many symmetric indefinite matrices. It may be desirable to restrict the use of this class to positive-definite matrices to guarantee accuracy.
|
||||||
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 FloatSymmetricLevinson.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:
IReadOnlyList<float32>
-
The right-hand side of the system.
Returns: FloatVector
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:
float32[]
-
The right-hand side of the system.
Returns: FloatVector
The solution vector.
|
Solve a symmetric square Toeplitz system with a right-side array. 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:
IROMatrix<float32>
-
The right-hand side of the system.
Returns: FloatMatrix
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 FloatSymmetricLevinson.IsSingular property be checked to see if the decomposition was completed, before attempting to obtain the upper triangle matrix.
|
Static members
| Static member | Description | ||||||
Full Usage:
FloatSymmetricLevinson.Inverse(T)
Parameters:
IReadOnlyList<float32>
-
The left-most column of the symmetric Toeplitz matrix.
Returns: FloatMatrix
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:
FloatSymmetricLevinson.Solve(T, Y)
Parameters:
IReadOnlyList<float32>
-
The left-most column of the Toeplitz matrix.
Y : IReadOnlyList<float32>
-
The right-side vector of the system.
Returns: FloatVector
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:
FloatSymmetricLevinson.Solve(T, Y)
Parameters:
IReadOnlyList<float32>
-
The left-most column of the Toeplitz matrix.
Y : IROMatrix<float32>
-
The right-side matrix of the system.
Returns: FloatMatrix
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:
FloatSymmetricLevinson.YuleWalker(R)
Parameters:
IReadOnlyList<float32>
-
The left-most column of the Toeplitz matrix.
Returns: FloatVector
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).
|