From 23810ef80e19e0319ebc97c1b35c41b5d23cd5d2 Mon Sep 17 00:00:00 2001 From: "Documenter.jl" Date: Thu, 4 Jan 2024 10:48:21 +0000 Subject: [PATCH] build based on d0e9cf5 --- dev/.documenter-siteinfo.json | 2 +- dev/condest.html | 8 ++++---- dev/index.html | 2 +- dev/iterative.html | 2 +- dev/lapackutil.html | 2 +- dev/lyapunov.html | 8 ++++---- dev/makeindex.html | 2 +- dev/meoperators.html | 8 ++++---- dev/meutil.html | 2 +- dev/plyapunov.html | 8 ++++---- dev/riccati.html | 10 +++++----- dev/sylvester.html | 14 +++++++------- dev/sylvkr.html | 6 +++--- 13 files changed, 37 insertions(+), 37 deletions(-) diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index 48419dc..13728ff 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.0","generation_timestamp":"2024-01-04T10:18:50","documenter_version":"1.2.1"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.0","generation_timestamp":"2024-01-04T10:48:16","documenter_version":"1.2.1"}} \ No newline at end of file diff --git a/dev/condest.html b/dev/condest.html index adf1460..87cb79e 100644 --- a/dev/condest.html +++ b/dev/condest.html @@ -9,7 +9,7 @@ 30.0 julia> opnorm1(invlyapop(A)) -3.7666666666666706source
MatrixEquations.opnorm1est — Function
γ = opnorm1est(op)

Compute γ, a lower bound of the 1-norm of the square linear operator op, using reverse communication based computations to evaluate op * x and op' * x. It is expected that in most cases $\gamma > \|op\|_1/10$, which is usually acceptable for estimating the condition numbers of linear operators.

Examples

julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
+3.7666666666666706
source
MatrixEquations.opnorm1est — Function
γ = opnorm1est(op)

Compute γ, a lower bound of the 1-norm of the square linear operator op, using reverse communication based computations to evaluate op * x and op' * x. It is expected that in most cases $\gamma > \|op\|_1/10$, which is usually acceptable for estimating the condition numbers of linear operators.

Examples

julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
 3×3 Array{Float64,2}:
  -6.0  -2.0   1.0
   5.0   1.0  -1.0
@@ -19,7 +19,7 @@
 18.0
 
 julia> opnorm1est(invlyapop(A))
-3.76666666666667
source
MatrixEquations.oprcondest — Function
rcond = oprcondest(op, opinv; exact = false)

Compute rcond, an estimation of the 1-norm reciprocal condition number of a linear operator op, where opinv is the inverse operator inv(op). The estimate is computed as $\text{rcond} = 1 / (\|op\|_1\|opinv\|_1)$, using estimates of the 1-norm, if exact = false, or computed exact values of the 1-norm, if exact = true. The exact = true option is not recommended for large order operators.

Note: No check is performed to verify that opinv = inv(op).

Examples

julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
+3.76666666666667
source
MatrixEquations.oprcondest — Function
rcond = oprcondest(op, opinv; exact = false)

Compute rcond, an estimation of the 1-norm reciprocal condition number of a linear operator op, where opinv is the inverse operator inv(op). The estimate is computed as $\text{rcond} = 1 / (\|op\|_1\|opinv\|_1)$, using estimates of the 1-norm, if exact = false, or computed exact values of the 1-norm, if exact = true. The exact = true option is not recommended for large order operators.

Note: No check is performed to verify that opinv = inv(op).

Examples

julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
 3×3 Array{Float64,2}:
  -6.0  -2.0   1.0
   5.0   1.0  -1.0
@@ -35,7 +35,7 @@
 0.008849557522123885
 
 julia> 1/opnorm1(lyapop(A))/opnorm1(invlyapop(A))
-0.008849557522123885
source
rcond = oprcondest(op; exact = false)

Compute rcond, an estimation of the 1-norm reciprocal condition number of a linear operator op, where op is one of the defined Lyapunov or Sylvester operators. The estimate is computed as $\text{rcond} = 1 / (\|op\|_1\|inv(op)\|_1)$, using estimates of the 1-norm, if exact = false, or computed exact values of the 1-norm, if exact = true. The exact = true option is not recommended for large order operators.

source
MatrixEquations.opsepest — Function
sep = opsepest(opinv; exact = false)

Compute sep, an estimation of the 1-norm separation of a linear operator op, where opinv is the inverse operator inv(op). The estimate is computed as $\text{sep} = 1 / \|opinv\|_1$ , using an estimate of the 1-norm, if exact = false, or the computed exact value of the 1-norm, if exact = true. The exact = true option is not recommended for large order operators.

The separation of the operator op is defined as

\[\text{sep} = \displaystyle\min_{X\neq 0} \frac{\|op*X\|}{\|X\|}.\]

An estimate of the reciprocal condition number of op can be computed as $\text{sep}/\|op\|_1$.

Example

julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
+0.008849557522123885
source
rcond = oprcondest(op; exact = false)

Compute rcond, an estimation of the 1-norm reciprocal condition number of a linear operator op, where op is one of the defined Lyapunov or Sylvester operators. The estimate is computed as $\text{rcond} = 1 / (\|op\|_1\|inv(op)\|_1)$, using estimates of the 1-norm, if exact = false, or computed exact values of the 1-norm, if exact = true. The exact = true option is not recommended for large order operators.

source
MatrixEquations.opsepest — Function
sep = opsepest(opinv; exact = false)

Compute sep, an estimation of the 1-norm separation of a linear operator op, where opinv is the inverse operator inv(op). The estimate is computed as $\text{sep} = 1 / \|opinv\|_1$ , using an estimate of the 1-norm, if exact = false, or the computed exact value of the 1-norm, if exact = true. The exact = true option is not recommended for large order operators.

The separation of the operator op is defined as

\[\text{sep} = \displaystyle\min_{X\neq 0} \frac{\|op*X\|}{\|X\|}.\]

An estimate of the reciprocal condition number of op can be computed as $\text{sep}/\|op\|_1$.

Example

julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
 3×3 Array{Float64,2}:
  -6.0  -2.0   1.0
   5.0   1.0  -1.0
@@ -51,4 +51,4 @@
 0.26548672566371656
 
 julia> 1/opnorm1(invlyapop(A))
-0.26548672566371656
source
+0.26548672566371656source diff --git a/dev/index.html b/dev/index.html index 40faf36..c46be34 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,2 +1,2 @@ -Overview · MatrixEquations.jl
GitHub

MatrixEquations.jl

DocBuild Code on Github.

This collection of Julia functions is an attemp to implement high performance numerical software to solve classes of Lyapunov, Sylvester and Riccati matrix equations at a performance level comparable with efficient structure exploiting Fortran implementations, as those available in the Systems and Control Library SLICOT. This goal has been fully achieved for Lyapunov and Sylvester equation solvers, for which the codes for both real and complex data perform at practically same performance level as similar functions available in the MATLAB Control System Toolbox (which rely on SLICOT).

The available functions in the MatrixEquations.jl package cover both standard and generalized continuous and discrete Lyapunov, Sylvester and Riccati equations for both real and complex data. The functions for the solution of Lyapunov and Sylvester equations rely on efficient structure exploiting solvers for which the input data are in Schur or generalized Schur forms. A comprehensive set of Lyapunov and Sylvester operators has been implemented, which allow the estimation of condition numbers of these operators and the iterative solution of various Lyapunov and Sylvester matrix equations using the conjugate gradient method. The implementation of Riccati equation solvers employ orthogonal Schur vectors based methods and their extensions to linear matrix pencil based reduction approaches. The calls of all functions with adjoint (in complex case) or transposed (in real case) arguments are fully supported by appropriate computational algorithms, thus the matrix copying operations are mostly avoided.

The current version of the package includes the following functions:

Solution of Lyapunov equations

FunctionDescription
lyapcSolution of the continuous Lyapunov equations
lyapdSolution of the discrete Lyapunov equations
tlyapcSolution of the continuous T-Lyapunov equations
hlyapcSolution of the continuous H-Lyapunov equations
tulyapc!Computation of the upper triangular solution of the continuous T-Lyapunov equation
hulyapc!Computation of the upper triangular solution of the continuous H-Lyapunov equation
plyapcSolution of the positive continuous Lyapunov equations
plyapdSolution of the positive discrete Lyapunov equations

Solution of algebraic Riccati equations

FunctionDescription
arecSolution of the continuous Riccati equations
garecSolution of the generalized continuous Riccati equation
aredSolution of the discrete Riccati equation
garedSolution of the generalized discrete Riccati equation

Solution of Sylvester equations and systems

FunctionDescription
sylvcSolution of the (continuous) Sylvester equations
sylvdSolution of the (discrete) Sylvester equations
gsylvSolution of the generalized Sylvester equations
sylvsysSolution of the Sylvester system of matrix equations
dsylvsysSolution of the dual Sylvester system of matrix equations

Iterative solution of linear matrix equations

FunctionDescription
lyapciIterative solution of continuous Lyapunov equations
lyapdiIterative solution of discrete Lyapunov equations
tlyapciIterative solution of the continuous T-Lyapunov equations
hlyapciIterative solution of the continuous H-Lyapunov equations
tulyapciIterative solution of the continuous T-Lyapunov equations with upper triangular solutions
hulyapciIterative solution of the continuous H-Lyapunov equations with upper triangular solutions
sylvciIterative solution of the (continuous) Sylvester equations
sylvdiIterative solution of the (discrete) Sylvester equations
gsylviIterative solution of the generalized Sylvester equations
gtsylviIterative solution of the generalized T-Sylvester equations
ghsylviIterative solution of the generalized H-Sylvester equations
cglsThe conjugate gradient method for nonsymmetric linear equations and least squares problems

Norm, condition and separation estimation of linear operators

FunctionDescription
opnorm1Computation of the 1-norm of a linear operator
opnorm1estEstimation of the 1-norm of a linear operator
oprcondestEstimation of the reciprocal 1-norm condition number of an operator
opsepestEstimation of the separation of a linear operator

The general solvers of Lyapunov and Sylvester equations rely on a set of specialized solvers for real or complex matrices in appropriate Schur forms. For testing purposes, a set of solvers for Sylvester equations has been implemented, which employ the Kronecker-product expansion of the equations. These solvers are not recommended for large order matrices. Based on the conjugate gradient method to solve linear systems or least-squares problems, several iterative solvers have been implemented, which are potentially applicable to solve linear matrix equations with large order dense and sparse matrices. The norms, reciprocal condition numbers and separations can be estimated for a comprehensive set of predefined Lyapunov and Sylvester operators. A complete list of implemented functions is available here.

Future plans

The collection of tools can be extended by adding new functionality, such as expert solvers, which additionally compute error bounds and condition estimates, or solvers for new classes of Sylvester-like equations or Riccati equations (as those arising in game-theoretic optimization problems). Further performance improvements are still possible by employing blocking based variants of solvers for Lyapunov and Sylvester equations.

Release Notes

Main developer

Andreas Varga

License: MIT (expat)

+Overview · MatrixEquations.jl
GitHub

MatrixEquations.jl

DocBuild Code on Github.

This collection of Julia functions is an attemp to implement high performance numerical software to solve classes of Lyapunov, Sylvester and Riccati matrix equations at a performance level comparable with efficient structure exploiting Fortran implementations, as those available in the Systems and Control Library SLICOT. This goal has been fully achieved for Lyapunov and Sylvester equation solvers, for which the codes for both real and complex data perform at practically same performance level as similar functions available in the MATLAB Control System Toolbox (which rely on SLICOT).

The available functions in the MatrixEquations.jl package cover both standard and generalized continuous and discrete Lyapunov, Sylvester and Riccati equations for both real and complex data. The functions for the solution of Lyapunov and Sylvester equations rely on efficient structure exploiting solvers for which the input data are in Schur or generalized Schur forms. A comprehensive set of Lyapunov and Sylvester operators has been implemented, which allow the estimation of condition numbers of these operators and the iterative solution of various Lyapunov and Sylvester matrix equations using the conjugate gradient method. The implementation of Riccati equation solvers employ orthogonal Schur vectors based methods and their extensions to linear matrix pencil based reduction approaches. The calls of all functions with adjoint (in complex case) or transposed (in real case) arguments are fully supported by appropriate computational algorithms, thus the matrix copying operations are mostly avoided.

The current version of the package includes the following functions:

Solution of Lyapunov equations

FunctionDescription
lyapcSolution of the continuous Lyapunov equations
lyapdSolution of the discrete Lyapunov equations
tlyapcSolution of the continuous T-Lyapunov equations
hlyapcSolution of the continuous H-Lyapunov equations
tulyapc!Computation of the upper triangular solution of the continuous T-Lyapunov equation
hulyapc!Computation of the upper triangular solution of the continuous H-Lyapunov equation
plyapcSolution of the positive continuous Lyapunov equations
plyapdSolution of the positive discrete Lyapunov equations

Solution of algebraic Riccati equations

FunctionDescription
arecSolution of the continuous Riccati equations
garecSolution of the generalized continuous Riccati equation
aredSolution of the discrete Riccati equation
garedSolution of the generalized discrete Riccati equation

Solution of Sylvester equations and systems

FunctionDescription
sylvcSolution of the (continuous) Sylvester equations
sylvdSolution of the (discrete) Sylvester equations
gsylvSolution of the generalized Sylvester equations
sylvsysSolution of the Sylvester system of matrix equations
dsylvsysSolution of the dual Sylvester system of matrix equations

Iterative solution of linear matrix equations

FunctionDescription
lyapciIterative solution of continuous Lyapunov equations
lyapdiIterative solution of discrete Lyapunov equations
tlyapciIterative solution of the continuous T-Lyapunov equations
hlyapciIterative solution of the continuous H-Lyapunov equations
tulyapciIterative solution of the continuous T-Lyapunov equations with upper triangular solutions
hulyapciIterative solution of the continuous H-Lyapunov equations with upper triangular solutions
sylvciIterative solution of the (continuous) Sylvester equations
sylvdiIterative solution of the (discrete) Sylvester equations
gsylviIterative solution of the generalized Sylvester equations
gtsylviIterative solution of the generalized T-Sylvester equations
ghsylviIterative solution of the generalized H-Sylvester equations
cglsThe conjugate gradient method for nonsymmetric linear equations and least squares problems

Norm, condition and separation estimation of linear operators

FunctionDescription
opnorm1Computation of the 1-norm of a linear operator
opnorm1estEstimation of the 1-norm of a linear operator
oprcondestEstimation of the reciprocal 1-norm condition number of an operator
opsepestEstimation of the separation of a linear operator

The general solvers of Lyapunov and Sylvester equations rely on a set of specialized solvers for real or complex matrices in appropriate Schur forms. For testing purposes, a set of solvers for Sylvester equations has been implemented, which employ the Kronecker-product expansion of the equations. These solvers are not recommended for large order matrices. Based on the conjugate gradient method to solve linear systems or least-squares problems, several iterative solvers have been implemented, which are potentially applicable to solve linear matrix equations with large order dense and sparse matrices. The norms, reciprocal condition numbers and separations can be estimated for a comprehensive set of predefined Lyapunov and Sylvester operators. A complete list of implemented functions is available here.

Future plans

The collection of tools can be extended by adding new functionality, such as expert solvers, which additionally compute error bounds and condition estimates, or solvers for new classes of Sylvester-like equations or Riccati equations (as those arising in game-theoretic optimization problems). Further performance improvements are still possible by employing blocking based variants of solvers for Lyapunov and Sylvester equations.

Release Notes

Main developer

Andreas Varga

License: MIT (expat)

diff --git a/dev/iterative.html b/dev/iterative.html index 6cea527..50dba6b 100644 --- a/dev/iterative.html +++ b/dev/iterative.html @@ -8,4 +8,4 @@ `info.resNE` - the relative residual for the normal equations `norm(A'*b - (A'*A + s*I)*x)/norm(A'*b)`; `info.iter` - the iteration number at which `x` was computed.

This function is a translation of the MATLAB implementation of CGLS, the conjugate gradient method for nonsymmetric linear equations and least squares problems https://web.stanford.edu/group/SOL/software/cgls/. The author of the code is Michael Saunders, with contributions from Per Christian Hansen, Folkert Bleichrodt and Christopher Fougner.

Note: Two alternative solvers lsqr and lsmr, available in the IterativeSolvers package, can also be employed. For example, the following call to lsqr can be alternatively used:

  using IterativeSolvers
-  lsqr(A, b; kwargs...) -> x[, history]

where kwargs contains solver-specific keyword arguments. A similar call to lsmr can be used.

source

Continuous-time Lyapunov and Lyapunov-like Matrix Equations

MatrixEquations.lyapci — Function
lyapci(A, C; abstol, reltol, maxiter) -> (X,info)

Compute for a square A and a hermitian/symmetric C a solution X of the continuous Lyapunov matrix equation

            A*X + X*A' + C = 0.

A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
lyapci(A, E, C; abstol, reltol, maxiter) -> (X,info)

Compute X, the solution of the generalized continuous Lyapunov equation

AXE' + EXA' + C = 0,

where A and E are square real or complex matrices and C is a square matrix. A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
MatrixEquations.lyapdi — Function
lyapdi(A, C; abstol, reltol, maxiter) -> (X,info)

Compute for a square A and a hermitian/symmetric C a solution X of the discrete Lyapunov matrix equation

            AXA' - X + C = 0.

A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
lyapdi(A, E, C; abstol, reltol, maxiter) -> (X,info)

Compute X, the solution of the generalized discrete Lyapunov equation

AXA' - EXE' + C = 0,

where A and E are square real or complex matrices and C is a square matrix. A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
MatrixEquations.tlyapci — Function
tlyapci(A, C, isig = +1; adj = false, abstol, reltol, maxiter) -> (X,info)

Compute a solution X of the continuous T-Lyapunov matrix equation

            A*X +isig*transpose(X)*transpose(A) = C   if adj = false,

or

            A*transpose(X) + isig*X*transpose(A) = C   if adj = true,

where for isig = 1, C is a symmetric matrix and for isig = -1, C is a skew-symmetric matrix.

For a matrix A, a least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined T-Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
MatrixEquations.hlyapci — Function
hlyapci(A, C, isig = +1; adj = false, abstol, reltol, maxiter) -> (X,info)

Compute a solution X of the continuous H-Lyapunov matrix equation

            A*X + isig*X'*A' = C   if adj = false,

or

            A*X' + isig*X*A' = C   if adj = true,

where for isig = 1, C is a hermitian matrix and for isig = -1, C is a skew-hermitian matrix.

For a matrix A, a least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined T-Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution. The keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
MatrixEquations.tulyapci — Function
tulyapci(U, Q; adj = false, abstol, reltol, maxiter) -> (X,info)

Compute for an upper triangular U and a symmetric Q an upper triangular solution X of the continuous T-Lyapunov matrix equation

  transpose(U)*X + transpose(X)*U = Q   if adj = false,

or 

  U*transpose(X) + X*transpose(U) = Q   if adj = true.

For a n×n upper triangular matrix U, a least-squares upper-triangular solution X is determined using a conjugate-gradient based iterative method applied to a suitably defined T-Lyapunov linear operator L:X -> Y, which maps upper triangular matrices X into upper triangular matrices Y, and the associated matrix M = Matrix(L) is $n(n+1)/2 \times n(n+1)/2$. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution. The keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
MatrixEquations.hulyapci — Function
hulyapci(U, Q; adj = false, abstol, reltol, maxiter) -> (X,info)

Compute for an upper triangular U and a hermitian Q an upper triangular solution X of the continuous H-Lyapunov matrix equation

            U'*X + X'*U = Q   if adj = false,

or

            U*X' + X*U' = Q    if adj = true.

For a n×n upper triangular matrix U, a least-squares upper-triangular solution X is determined using a conjugate-gradient based iterative method applied to a suitably defined T-Lyapunov linear operator L:X -> Y, which maps upper triangular matrices X into upper triangular matrices Y, and the associated matrix M = Matrix(L) is $n(n+1)/2 \times n(n+1)/2$. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution. The keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source

Sylvester and Sylvester-like Matrix Equations

MatrixEquations.sylvci — Function
X = sylvci(A,B,C)

Solve the continuous Sylvester matrix equation

            AX + XB = C ,

where A and B are square matrices.

A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
MatrixEquations.sylvdi — Function
X = sylvdi(A,B,C)

Solve the discrete Sylvester matrix equation

            AXB + X = C ,

where A and B are square matrices.

A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
MatrixEquations.gsylvi — Function
X = gsylvi(A,B,C,D,E)

Solve the generalized Sylvester matrix equation

AXB + CXD = E ,

where A, B, C and D are square matrices.

A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
MatrixEquations.gtsylvi — Function
 gtsylvi(A, B, C, D, E; mx, nx, abstol, reltol, maxiter) -> (X,info)

Compute a solution X of the generalized T-Sylvester matrix equation

  ∑ A_i*X*B_i + ∑ C_j*transpose(X)*D_j = E,

where A_i and C_j are matrices having the same row dimension equal to the row dimension of E and B_i and D_j are matrices having the same column dimension equal to the column dimension of E. A_i and B_i are contained in the k-vectors of matrices A and B, respectively, and C_j and D_j are contained in the l-vectors of matrices C and D, respectively. Any of the component matrices can be given as an UniformScaling. The keyword parameters mx and nx can be used to specify the row and column dimensions of X, if they cannot be inferred from the data.

A least-squares solution X is determined using a conjugate-gradient based iterative method applied to a suitably defined T-Sylvester linear operator L:X -> Y such that L(X) = E or norm(L(X) - E) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

Note: For the derivation of the adjoint equation see reference [1], which also served as motivation to implement a general linear matrix equation solver in Julia.

[1] Uhlig, F., Xu, A.B. Iterative optimal solutions of linear matrix equations for hyperspectral and multispectral image fusing. Calcolo 60, 26 (2023). https://doi.org/10.1007/s10092-023-00514-8

source
MatrixEquations.ghsylvi — Function
 ghsylvi(A, B, C, D, E; mx, nx, abstol, reltol, maxiter) -> (X,info)

Compute a solution X of the generalized H-Sylvester matrix equation

  ∑ A_i*X*B_i + ∑ C_j*X'*D_j = E,

where A_i and C_j are matrices having the same row dimension equal to the row dimension of E and B_i and D_j are matrices having the same column dimension equal to the column dimension of E. A_i and B_i are contained in the k-vectors of matrices A and B, respectively, and C_j and D_j are contained in the l-vectors of matrices C and D, respectively. Any of the component matrices can be given as an UniformScaling. The keyword parameters mx and nx can be used to specify the row and column dimensions of X, if they cannot be inferred from the data.

A least-squares solution X is determined using a conjugate-gradient based iterative method applied to a suitably defined T-Sylvester linear operator L:X -> Y such that L(X) = E or norm(L(X) - E) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

Note: For the derivation of the adjoint equation see reference [1], which also served as motivation to implement a general linear matrix equation solver in Julia.

[1] Uhlig, F., Xu, A.B. Iterative optimal solutions of linear matrix equations for hyperspectral and multispectral image fusing. Calcolo 60, 26 (2023). https://doi.org/10.1007/s10092-023-00514-8

source
+ lsqr(A, b; kwargs...) -> x[, history]

where kwargs contains solver-specific keyword arguments. A similar call to lsmr can be used.

source

Continuous-time Lyapunov and Lyapunov-like Matrix Equations

MatrixEquations.lyapci — Function
lyapci(A, C; abstol, reltol, maxiter) -> (X,info)

Compute for a square A and a hermitian/symmetric C a solution X of the continuous Lyapunov matrix equation

            A*X + X*A' + C = 0.

A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
lyapci(A, E, C; abstol, reltol, maxiter) -> (X,info)

Compute X, the solution of the generalized continuous Lyapunov equation

AXE' + EXA' + C = 0,

where A and E are square real or complex matrices and C is a square matrix. A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
MatrixEquations.lyapdi — Function
lyapdi(A, C; abstol, reltol, maxiter) -> (X,info)

Compute for a square A and a hermitian/symmetric C a solution X of the discrete Lyapunov matrix equation

            AXA' - X + C = 0.

A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
lyapdi(A, E, C; abstol, reltol, maxiter) -> (X,info)

Compute X, the solution of the generalized discrete Lyapunov equation

AXA' - EXE' + C = 0,

where A and E are square real or complex matrices and C is a square matrix. A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
MatrixEquations.tlyapci — Function
tlyapci(A, C, isig = +1; adj = false, abstol, reltol, maxiter) -> (X,info)

Compute a solution X of the continuous T-Lyapunov matrix equation

            A*X +isig*transpose(X)*transpose(A) = C   if adj = false,

or

            A*transpose(X) + isig*X*transpose(A) = C   if adj = true,

where for isig = 1, C is a symmetric matrix and for isig = -1, C is a skew-symmetric matrix.

For a matrix A, a least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined T-Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
MatrixEquations.hlyapci — Function
hlyapci(A, C, isig = +1; adj = false, abstol, reltol, maxiter) -> (X,info)

Compute a solution X of the continuous H-Lyapunov matrix equation

            A*X + isig*X'*A' = C   if adj = false,

or

            A*X' + isig*X*A' = C   if adj = true,

where for isig = 1, C is a hermitian matrix and for isig = -1, C is a skew-hermitian matrix.

For a matrix A, a least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined T-Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution. The keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
MatrixEquations.tulyapci — Function
tulyapci(U, Q; adj = false, abstol, reltol, maxiter) -> (X,info)

Compute for an upper triangular U and a symmetric Q an upper triangular solution X of the continuous T-Lyapunov matrix equation

  transpose(U)*X + transpose(X)*U = Q   if adj = false,

or 

  U*transpose(X) + X*transpose(U) = Q   if adj = true.

For a n×n upper triangular matrix U, a least-squares upper-triangular solution X is determined using a conjugate-gradient based iterative method applied to a suitably defined T-Lyapunov linear operator L:X -> Y, which maps upper triangular matrices X into upper triangular matrices Y, and the associated matrix M = Matrix(L) is $n(n+1)/2 \times n(n+1)/2$. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution. The keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
MatrixEquations.hulyapci — Function
hulyapci(U, Q; adj = false, abstol, reltol, maxiter) -> (X,info)

Compute for an upper triangular U and a hermitian Q an upper triangular solution X of the continuous H-Lyapunov matrix equation

            U'*X + X'*U = Q   if adj = false,

or

            U*X' + X*U' = Q    if adj = true.

For a n×n upper triangular matrix U, a least-squares upper-triangular solution X is determined using a conjugate-gradient based iterative method applied to a suitably defined T-Lyapunov linear operator L:X -> Y, which maps upper triangular matrices X into upper triangular matrices Y, and the associated matrix M = Matrix(L) is $n(n+1)/2 \times n(n+1)/2$. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution. The keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source

Sylvester and Sylvester-like Matrix Equations

MatrixEquations.sylvci — Function
X = sylvci(A,B,C)

Solve the continuous Sylvester matrix equation

            AX + XB = C ,

where A and B are square matrices.

A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
MatrixEquations.sylvdi — Function
X = sylvdi(A,B,C)

Solve the discrete Sylvester matrix equation

            AXB + X = C ,

where A and B are square matrices.

A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
MatrixEquations.gsylvi — Function
X = gsylvi(A,B,C,D,E)

Solve the generalized Sylvester matrix equation

AXB + CXD = E ,

where A, B, C and D are square matrices.

A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
MatrixEquations.gtsylvi — Function
 gtsylvi(A, B, C, D, E; mx, nx, abstol, reltol, maxiter) -> (X,info)

Compute a solution X of the generalized T-Sylvester matrix equation

  ∑ A_i*X*B_i + ∑ C_j*transpose(X)*D_j = E,

where A_i and C_j are matrices having the same row dimension equal to the row dimension of E and B_i and D_j are matrices having the same column dimension equal to the column dimension of E. A_i and B_i are contained in the k-vectors of matrices A and B, respectively, and C_j and D_j are contained in the l-vectors of matrices C and D, respectively. Any of the component matrices can be given as an UniformScaling. The keyword parameters mx and nx can be used to specify the row and column dimensions of X, if they cannot be inferred from the data.

A least-squares solution X is determined using a conjugate-gradient based iterative method applied to a suitably defined T-Sylvester linear operator L:X -> Y such that L(X) = E or norm(L(X) - E) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

Note: For the derivation of the adjoint equation see reference [1], which also served as motivation to implement a general linear matrix equation solver in Julia.

[1] Uhlig, F., Xu, A.B. Iterative optimal solutions of linear matrix equations for hyperspectral and multispectral image fusing. Calcolo 60, 26 (2023). https://doi.org/10.1007/s10092-023-00514-8

source
MatrixEquations.ghsylvi — Function
 ghsylvi(A, B, C, D, E; mx, nx, abstol, reltol, maxiter) -> (X,info)

Compute a solution X of the generalized H-Sylvester matrix equation

  ∑ A_i*X*B_i + ∑ C_j*X'*D_j = E,

where A_i and C_j are matrices having the same row dimension equal to the row dimension of E and B_i and D_j are matrices having the same column dimension equal to the column dimension of E. A_i and B_i are contained in the k-vectors of matrices A and B, respectively, and C_j and D_j are contained in the l-vectors of matrices C and D, respectively. Any of the component matrices can be given as an UniformScaling. The keyword parameters mx and nx can be used to specify the row and column dimensions of X, if they cannot be inferred from the data.

A least-squares solution X is determined using a conjugate-gradient based iterative method applied to a suitably defined T-Sylvester linear operator L:X -> Y such that L(X) = E or norm(L(X) - E) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

Note: For the derivation of the adjoint equation see reference [1], which also served as motivation to implement a general linear matrix equation solver in Julia.

[1] Uhlig, F., Xu, A.B. Iterative optimal solutions of linear matrix equations for hyperspectral and multispectral image fusing. Calcolo 60, 26 (2023). https://doi.org/10.1007/s10092-023-00514-8

source
diff --git a/dev/lapackutil.html b/dev/lapackutil.html index 0443772..83182c6 100644 --- a/dev/lapackutil.html +++ b/dev/lapackutil.html @@ -1,4 +1,4 @@ Lapack Utilities · MatrixEquations.jl
GitHub

Lapack Utilities

MatrixEquations.LapackUtil.tgsyl! — Function
tgsyl!(A, B, C, D, E, F) -> (C, F, scale)

Solve the Sylvester system of matrix equations

  AX - YB = scale*C
   DX - YE = scale*F ,

where X and Y are unknown matrices, the pairs (A, D), (B, E) and (C, F) have the same sizes, and the pairs (A, D) and (B, E) are in generalized (real) Schur canonical form, i.e. A, B are upper quasi triangular and D, E are upper triangular. Returns X (overwriting C), Y (overwriting F) and scale.

tgsyl!(trans, A, B, C, D, E, F) -> (C, F, scale)

Solve for trans = 'T' and real matrices or for trans = 'C' and complex matrices, the (adjoint) Sylvester system of matrix equations

  A'X + D'Y = scale*C
-  XB' + YE' = scale*(-F) .

tgsyl!('N', A, B, C, D, E, F) corresponds to the call tgsyl!(A, B, C, D, E, F).

Interface to the LAPACK subroutines DTGSYL/STGSYL/ZTGSYL/CTGSYL.

source
MatrixEquations.LapackUtil.lanv2 — Function
lanv2(A, B, C, D) -> (RT1R, RT1I, RT2R, RT2I, CS, SN)

Compute the Schur factorization of a real 2-by-2 nonsymmetric matrix [A,B;C,D] in standard form. A, B, C, D are overwritten on output by the corresponding elements of the standardised Schur form. RT1R+im*RT1I and RT2R+im*RT2I are the resulting eigenvalues. CS and SN are the parameters of the rotation matrix. Interface to the LAPACK subroutines DLANV2/SLANV2.

source
MatrixEquations.LapackUtil.lag2 — Function
lag2(A, B, SAFMIN) -> (SCALE1, SCALE2, WR1, WR2, WI)

Compute the eigenvalues of a 2-by-2 generalized real eigenvalue problem for the matrix pair (A,B), with scaling as necessary to avoid over-/underflow. SAFMIN is the smallest positive number s.t. 1/SAFMIN does not overflow. If WI = 0, WR1/SCALE1 and WR2/SCALE2 are the resulting real eigenvalues, while if WI <> 0, then (WR1+/-im*WI)/SCALE1 are the resulting complex eigenvalues. Interface to the LAPACK subroutines DLAG2/SLAG2.

source
MatrixEquations.LapackUtil.ladiv — Function
ladiv(A, B, C, D) -> (P, Q)

Perform the complex division in real arithmetic

$P + iQ = \displaystyle\frac{A+iB}{C+iD}$

by avoiding unnecessary overflow. Interface to the LAPACK subroutines DLADIV/SLADIV.

source
MatrixEquations.LapackUtil.lacn2! — Function
lacn2!(V, X, ISGN, EST, KASE, ISAVE ) -> (EST, KASE )

Estimate the 1-norm of a real linear operator A, using reverse communication by applying the operator or its transpose/adjoint to a vector. KASE is a parameter to control the norm evaluation process as follows. On the initial call, KASE should be 0. On an intermediate return, KASE will be 1 or 2, indicating whether the real vector X should be overwritten by A * X or A' * X at the next call. On the final return, KASE will again be 0 and EST is an estimate (a lower bound) for the 1-norm of A. V is a real work vector, ISGN is an integer work vector and ISAVE is a 3-dimensional integer vector used to save information between the calls. Interface to the LAPACK subroutines DLACN2/SLACN2.

source
lacn2!(V, X, EST, KASE, ISAVE ) -> (EST, KASE )

Estimate the 1-norm of a complex linear operator A, using reverse communication by applying the operator or its adjoint to a vector. KASE is a parameter to control the norm evaluation process as follows. On the initial call, KASE should be 0. On an intermediate return, KASE will be 1 or 2, indicating whether the complex vector X should be overwritten by A * X or A' * X at the next call. On the final return, KASE will again be 0 and EST is an estimate (a lower bound) for the 1-norm of A. V is a complex work vector and ISAVE is a 3-dimensional integer vector used to save information between the calls. Interface to the LAPACK subroutines ZLACN2/CLACN2.

source
+ XB' + YE' = scale*(-F) .

tgsyl!('N', A, B, C, D, E, F) corresponds to the call tgsyl!(A, B, C, D, E, F).

Interface to the LAPACK subroutines DTGSYL/STGSYL/ZTGSYL/CTGSYL.

source
MatrixEquations.LapackUtil.lanv2 — Function
lanv2(A, B, C, D) -> (RT1R, RT1I, RT2R, RT2I, CS, SN)

Compute the Schur factorization of a real 2-by-2 nonsymmetric matrix [A,B;C,D] in standard form. A, B, C, D are overwritten on output by the corresponding elements of the standardised Schur form. RT1R+im*RT1I and RT2R+im*RT2I are the resulting eigenvalues. CS and SN are the parameters of the rotation matrix. Interface to the LAPACK subroutines DLANV2/SLANV2.

source
MatrixEquations.LapackUtil.lag2 — Function
lag2(A, B, SAFMIN) -> (SCALE1, SCALE2, WR1, WR2, WI)

Compute the eigenvalues of a 2-by-2 generalized real eigenvalue problem for the matrix pair (A,B), with scaling as necessary to avoid over-/underflow. SAFMIN is the smallest positive number s.t. 1/SAFMIN does not overflow. If WI = 0, WR1/SCALE1 and WR2/SCALE2 are the resulting real eigenvalues, while if WI <> 0, then (WR1+/-im*WI)/SCALE1 are the resulting complex eigenvalues. Interface to the LAPACK subroutines DLAG2/SLAG2.

source
MatrixEquations.LapackUtil.ladiv — Function
ladiv(A, B, C, D) -> (P, Q)

Perform the complex division in real arithmetic

$P + iQ = \displaystyle\frac{A+iB}{C+iD}$

by avoiding unnecessary overflow. Interface to the LAPACK subroutines DLADIV/SLADIV.

source
MatrixEquations.LapackUtil.lacn2! — Function
lacn2!(V, X, ISGN, EST, KASE, ISAVE ) -> (EST, KASE )

Estimate the 1-norm of a real linear operator A, using reverse communication by applying the operator or its transpose/adjoint to a vector. KASE is a parameter to control the norm evaluation process as follows. On the initial call, KASE should be 0. On an intermediate return, KASE will be 1 or 2, indicating whether the real vector X should be overwritten by A * X or A' * X at the next call. On the final return, KASE will again be 0 and EST is an estimate (a lower bound) for the 1-norm of A. V is a real work vector, ISGN is an integer work vector and ISAVE is a 3-dimensional integer vector used to save information between the calls. Interface to the LAPACK subroutines DLACN2/SLACN2.

source
lacn2!(V, X, EST, KASE, ISAVE ) -> (EST, KASE )

Estimate the 1-norm of a complex linear operator A, using reverse communication by applying the operator or its adjoint to a vector. KASE is a parameter to control the norm evaluation process as follows. On the initial call, KASE should be 0. On an intermediate return, KASE will be 1 or 2, indicating whether the complex vector X should be overwritten by A * X or A' * X at the next call. On the final return, KASE will again be 0 and EST is an estimate (a lower bound) for the 1-norm of A. V is a complex work vector and ISAVE is a 3-dimensional integer vector used to save information between the calls. Interface to the LAPACK subroutines ZLACN2/CLACN2.

source
diff --git a/dev/lyapunov.html b/dev/lyapunov.html index 9b269d7..aab7159 100644 --- a/dev/lyapunov.html +++ b/dev/lyapunov.html @@ -17,7 +17,7 @@ julia> A*X + X*A' + C 2×2 Array{Float64,2}: -8.88178e-16 2.22045e-16 - 2.22045e-16 -4.44089e-16source
X = lyapc(A, E, C)

Compute X, the solution of the generalized continuous Lyapunov equation

 AXE' + EXA' + C = 0,

where A and E are square real or complex matrices and C is a square matrix. The pencil A-λE must not have two eigenvalues α and β such that α+β = 0. The solution X is symmetric or hermitian if C is a symmetric or hermitian.

The following particular cases are also adressed:

X = lyapc(A,β*I,C)  or  X = lyapc(A,β,C)

Solve the matrix equation AXβ' + βXA' + C = 0.

X = lyapc(α*I,E,C)  or  X = lyapc(α,E,C)

Solve the matrix equation αXE' + EXα' + C = 0.

X = lyapc(α*I,β*I,C)  or  X = lyapc(α,β,C)

Solve the matrix equation (αβ'+α'β)X + C = 0.

x = lyapc(α,β,γ)

Solve the equation (αβ'+α'β)x + γ = 0.

Example

julia> A = [3. 4.; 5. 6.]
+  2.22045e-16  -4.44089e-16
source
X = lyapc(A, E, C)

Compute X, the solution of the generalized continuous Lyapunov equation

 AXE' + EXA' + C = 0,

where A and E are square real or complex matrices and C is a square matrix. The pencil A-λE must not have two eigenvalues α and β such that α+β = 0. The solution X is symmetric or hermitian if C is a symmetric or hermitian.

The following particular cases are also adressed:

X = lyapc(A,β*I,C)  or  X = lyapc(A,β,C)

Solve the matrix equation AXβ' + βXA' + C = 0.

X = lyapc(α*I,E,C)  or  X = lyapc(α,E,C)

Solve the matrix equation αXE' + EXα' + C = 0.

X = lyapc(α*I,β*I,C)  or  X = lyapc(α,β,C)

Solve the matrix equation (αβ'+α'β)X + C = 0.

x = lyapc(α,β,γ)

Solve the equation (αβ'+α'β)x + γ = 0.

Example

julia> A = [3. 4.; 5. 6.]
 2×2 Array{Float64,2}:
  3.0  4.0
  5.0  6.0
@@ -40,7 +40,7 @@
 julia> A*X*E' + E*X*A' + C
 2×2 Array{Float64,2}:
  -5.32907e-15  -2.66454e-15
- -4.44089e-15   0.0
source
MatrixEquations.lyapcs! — Function
lyapcs!(A,C;adj = false)

Solve the continuous Lyapunov matrix equation

            op(A)X + Xop(A)' + C = 0,

where op(A) = A if adj = false and op(A) = A' if adj = true. A is a square real matrix in a real Schur form, or a square complex matrix in a complex Schur form and C is a symmetric or hermitian matrix. A must not have two eigenvalues α and β such that α+β = 0. C contains on output the solution X.

source
lyapcs!(A, E, C; adj = false)

Solve the generalized continuous Lyapunov matrix equation

            op(A)Xop(E)' + op(E)Xop(A)' + C = 0,

where op(A) = A and op(E) = E if adj = false and op(A) = A' and op(E) = E' if adj = true. The pair (A,E) is in a generalized real or complex Schur form and C is a symmetric or hermitian matrix. The pencil A-λE must not have two eigenvalues α and β such that α+β = 0. The computed symmetric or hermitian solution X is contained in C.

source
MatrixEquations.tlyapc — Function
X = tlyapc(A,C, isig = 1; fast = true, atol::Real=0, rtol::Real=atol>0 ? 0 : N*ϵ)

Compute for isig = ±1 a solution of the continuous T-Lyapunov matrix equation

            A*X + isig*transpose(X)*transpose(A) + C = 0

using explicit formulas based on full rank factorizations of A (see [1] and [2]). A and C are m×n and m×m matrices, respectively, and X is an n×m matrix. The matrix C must be symmetric if isig = 1 and skew-symmetric if isig = -1. atol and rtol are the absolute and relative tolerances, respectively, used for rank computation. The default relative tolerance is N*ϵ, where N = min(m,n) and ϵ is the machine precision of the element type of A.

The underlying rank revealing factorization of A is the QR-factorization with column pivoting, if fast = true (default), or the more reliable SVD-decomposition, if fast = false.

[1] H. W. Braden. The equations A^T X ± X^T A = B. SIAM J. Matrix Anal. Appl., 20(2):295–302, 1998.

[2] C.-Y. Chiang, E. K.-W. Chu, W.-W. Lin, On the ★-Sylvester equation AX ± X^★ B^★ = C. Applied Mathematics and Computation, 218:8393–8407, 2012.

source
MatrixEquations.hlyapc — Function
X = hlyapc(A,C, isig = 1; atol::Real=0, rtol::Real=atol>0 ? 0 : N*ϵ)

Compute for isig = ±1 a solution of the continuous H-Lyapunov matrix equation

            A*X + isig*adjoint(X)*adjoint(A) + C = 0

using explicit formulas based on full rank factorizations of A (see [1] and [2]). A and C are m×n and m×m matrices, respectively, and X is an n×m matrix. The matrix C must be hermitian if isig = 1 and skew-hermitian if isig = -1. atol and rtol are the absolute and relative tolerances, respectively, used for rank computation. The default relative tolerance is N*ϵ, where N = min(m,n) and ϵ is the machine precision of the element type of A.

The underlying rank revealing factorization of Ae (or of A if real) is the QR-factorization with column pivoting, if fast = true (default), or the more reliable SVD-decomposition, if fast = false.

[1] H. W. Braden. The equations A^T X ± X^T A = B. SIAM J. Matrix Anal. Appl., 20(2):295–302, 1998.

[2] C.-Y. Chiang, E. K.-W. Chu, W.-W. Lin, On the ★-Sylvester equation AX ± X^★ B^★ = C. Applied Mathematics and Computation, 218:8393–8407, 2012.

source
MatrixEquations.tulyapc! — Function
X = tulyapc!(U, Q; adj = false)

Compute for an upper triangular U and a symmetric Q an upper triangular solution X of the continuous T-Lyapunov matrix equation

            transpose(U)*X + transpose(X)*U = Q  if adj = false,

or

            U*transpose(X) + X*transpose(U) = Q   if adj = true.

The solution X overwrites the matrix Q, while U is unchanged.

If U is nonsingular, a backward elimination method is used, if adj = false, or a forward elimination method is used if adj = true, by adapting the approach employed in the function dU_from_dQ! of the DiffOpt.jl package. For a n×n singular U, a least-squares upper-triangular solution X is determined using a conjugate-gradient based iterative method applied to a suitably defined T-Lyapunov linear operator L:X -> Y, which maps upper triangular matrices X into upper triangular matrices Y, and the associated matrix M = Matrix(L) is $n(n+1)/2 \times n(n+1)/2$. In this case, the keyword argument tol (default: tol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
MatrixEquations.hulyapc! — Function
X = hulyapc!(U, Q; adj = false)

Compute for an upper triangular U and a hermitian Q an upper triangular solution X of the continuous H-Lyapunov matrix equation

            U'*X + X'*U = Q   if adj = false,

or

            U*X' + X*U' = Q   if adj = true.

The solution X overwrites the matrix Q, while U is unchanged.

If U is nonsingular, a backward elimination method is used, if adj = false, or a forward elimination method is used if adj = true, by adapting the approach employed in the function dU_from_dQ! of the DiffOpt.jl package. For a n×n singular U, a least-squares upper-triangular solution X is determined using a conjugate-gradient based iterative method applied to a suitably defined T-Lyapunov linear operator L:X -> Y, which maps upper triangular matrices X into upper triangular matrices Y, and the associated matrix M = Matrix(L) is $n(n+1)/2 \times n(n+1)/2$. In this case, the keyword argument tol (default: tol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source

Discrete-time Lyapunov (Stein) Matrix Equations

MatrixEquations.lyapd — Function
X = lyapd(A, C)

Compute X, the solution of the discrete Lyapunov equation

   AXA' - X + C = 0,

where A is a square real or complex matrix and C is a square matrix. A must not have two eigenvalues α and β such that αβ = 1. The solution X is symmetric or hermitian if C is a symmetric or hermitian.

The following particular cases are also adressed:

X = lyapd(α*I,C)  or  X = lyapd(α,C)

Solve the matrix equation (αα'-1)X + C = 0.

x = lyapd(α,γ)

Solve the equation (αα'-1)x + γ = 0.

Example

julia> A = [3. 4.; 5. 6.]
+ -4.44089e-15   0.0
source
MatrixEquations.lyapcs! — Function
lyapcs!(A,C;adj = false)

Solve the continuous Lyapunov matrix equation

            op(A)X + Xop(A)' + C = 0,

where op(A) = A if adj = false and op(A) = A' if adj = true. A is a square real matrix in a real Schur form, or a square complex matrix in a complex Schur form and C is a symmetric or hermitian matrix. A must not have two eigenvalues α and β such that α+β = 0. C contains on output the solution X.

source
lyapcs!(A, E, C; adj = false)

Solve the generalized continuous Lyapunov matrix equation

            op(A)Xop(E)' + op(E)Xop(A)' + C = 0,

where op(A) = A and op(E) = E if adj = false and op(A) = A' and op(E) = E' if adj = true. The pair (A,E) is in a generalized real or complex Schur form and C is a symmetric or hermitian matrix. The pencil A-λE must not have two eigenvalues α and β such that α+β = 0. The computed symmetric or hermitian solution X is contained in C.

source
MatrixEquations.tlyapc — Function
X = tlyapc(A,C, isig = 1; fast = true, atol::Real=0, rtol::Real=atol>0 ? 0 : N*ϵ)

Compute for isig = ±1 a solution of the continuous T-Lyapunov matrix equation

            A*X + isig*transpose(X)*transpose(A) + C = 0

using explicit formulas based on full rank factorizations of A (see [1] and [2]). A and C are m×n and m×m matrices, respectively, and X is an n×m matrix. The matrix C must be symmetric if isig = 1 and skew-symmetric if isig = -1. atol and rtol are the absolute and relative tolerances, respectively, used for rank computation. The default relative tolerance is N*ϵ, where N = min(m,n) and ϵ is the machine precision of the element type of A.

The underlying rank revealing factorization of A is the QR-factorization with column pivoting, if fast = true (default), or the more reliable SVD-decomposition, if fast = false.

[1] H. W. Braden. The equations A^T X ± X^T A = B. SIAM J. Matrix Anal. Appl., 20(2):295–302, 1998.

[2] C.-Y. Chiang, E. K.-W. Chu, W.-W. Lin, On the ★-Sylvester equation AX ± X^★ B^★ = C. Applied Mathematics and Computation, 218:8393–8407, 2012.

source
MatrixEquations.hlyapc — Function
X = hlyapc(A,C, isig = 1; atol::Real=0, rtol::Real=atol>0 ? 0 : N*ϵ)

Compute for isig = ±1 a solution of the continuous H-Lyapunov matrix equation

            A*X + isig*adjoint(X)*adjoint(A) + C = 0

using explicit formulas based on full rank factorizations of A (see [1] and [2]). A and C are m×n and m×m matrices, respectively, and X is an n×m matrix. The matrix C must be hermitian if isig = 1 and skew-hermitian if isig = -1. atol and rtol are the absolute and relative tolerances, respectively, used for rank computation. The default relative tolerance is N*ϵ, where N = min(m,n) and ϵ is the machine precision of the element type of A.

The underlying rank revealing factorization of Ae (or of A if real) is the QR-factorization with column pivoting, if fast = true (default), or the more reliable SVD-decomposition, if fast = false.

[1] H. W. Braden. The equations A^T X ± X^T A = B. SIAM J. Matrix Anal. Appl., 20(2):295–302, 1998.

[2] C.-Y. Chiang, E. K.-W. Chu, W.-W. Lin, On the ★-Sylvester equation AX ± X^★ B^★ = C. Applied Mathematics and Computation, 218:8393–8407, 2012.

source
MatrixEquations.tulyapc! — Function
X = tulyapc!(U, Q; adj = false)

Compute for an upper triangular U and a symmetric Q an upper triangular solution X of the continuous T-Lyapunov matrix equation

            transpose(U)*X + transpose(X)*U = Q  if adj = false,

or

            U*transpose(X) + X*transpose(U) = Q   if adj = true.

The solution X overwrites the matrix Q, while U is unchanged.

If U is nonsingular, a backward elimination method is used, if adj = false, or a forward elimination method is used if adj = true, by adapting the approach employed in the function dU_from_dQ! of the DiffOpt.jl package. For a n×n singular U, a least-squares upper-triangular solution X is determined using a conjugate-gradient based iterative method applied to a suitably defined T-Lyapunov linear operator L:X -> Y, which maps upper triangular matrices X into upper triangular matrices Y, and the associated matrix M = Matrix(L) is $n(n+1)/2 \times n(n+1)/2$. In this case, the keyword argument tol (default: tol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source
MatrixEquations.hulyapc! — Function
X = hulyapc!(U, Q; adj = false)

Compute for an upper triangular U and a hermitian Q an upper triangular solution X of the continuous H-Lyapunov matrix equation

            U'*X + X'*U = Q   if adj = false,

or

            U*X' + X*U' = Q   if adj = true.

The solution X overwrites the matrix Q, while U is unchanged.

If U is nonsingular, a backward elimination method is used, if adj = false, or a forward elimination method is used if adj = true, by adapting the approach employed in the function dU_from_dQ! of the DiffOpt.jl package. For a n×n singular U, a least-squares upper-triangular solution X is determined using a conjugate-gradient based iterative method applied to a suitably defined T-Lyapunov linear operator L:X -> Y, which maps upper triangular matrices X into upper triangular matrices Y, and the associated matrix M = Matrix(L) is $n(n+1)/2 \times n(n+1)/2$. In this case, the keyword argument tol (default: tol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

source

Discrete-time Lyapunov (Stein) Matrix Equations

MatrixEquations.lyapd — Function
X = lyapd(A, C)

Compute X, the solution of the discrete Lyapunov equation

   AXA' - X + C = 0,

where A is a square real or complex matrix and C is a square matrix. A must not have two eigenvalues α and β such that αβ = 1. The solution X is symmetric or hermitian if C is a symmetric or hermitian.

The following particular cases are also adressed:

X = lyapd(α*I,C)  or  X = lyapd(α,C)

Solve the matrix equation (αα'-1)X + C = 0.

x = lyapd(α,γ)

Solve the equation (αα'-1)x + γ = 0.

Example

julia> A = [3. 4.; 5. 6.]
 2×2 Array{Float64,2}:
  3.0  4.0
  5.0  6.0
@@ -58,7 +58,7 @@
 julia> A*X*A' - X + C
 2×2 Array{Float64,2}:
  5.55112e-16  6.66134e-16
- 2.22045e-16  4.44089e-16
source
X = lyapd(A, E, C)

Compute X, the solution of the generalized discrete Lyapunov equation

     AXA' - EXE' + C = 0,

where A and E are square real or complex matrices and C is a square matrix. The pencil A-λE must not have two eigenvalues α and β such that αβ = 1. The solution X is symmetric or hermitian if C is a symmetric or hermitian.

The following particular cases are also adressed:

X = lyapd(A,β*I,C)  or  X = lyapd(A,β,C)

Solve the matrix equation AXA' - βXβ' + C = 0.

X = lyapd(α*I,E,C)  or  X = lyapd(α,E,C)

Solve the matrix equation αXα' - EXE' + C = 0.

X = lyapd(α*I,β*I,C)  or  X = lyapd(α,β,C)

Solve the matrix equation (αα'-ββ')X + C = 0.

x = lyapd(α,β,γ)

Solve the equation (αα'-ββ')x + γ = 0.

Example

julia> A = [3. 4.; 5. 6.]
+ 2.22045e-16  4.44089e-16
source
X = lyapd(A, E, C)

Compute X, the solution of the generalized discrete Lyapunov equation

     AXA' - EXE' + C = 0,

where A and E are square real or complex matrices and C is a square matrix. The pencil A-λE must not have two eigenvalues α and β such that αβ = 1. The solution X is symmetric or hermitian if C is a symmetric or hermitian.

The following particular cases are also adressed:

X = lyapd(A,β*I,C)  or  X = lyapd(A,β,C)

Solve the matrix equation AXA' - βXβ' + C = 0.

X = lyapd(α*I,E,C)  or  X = lyapd(α,E,C)

Solve the matrix equation αXα' - EXE' + C = 0.

X = lyapd(α*I,β*I,C)  or  X = lyapd(α,β,C)

Solve the matrix equation (αα'-ββ')X + C = 0.

x = lyapd(α,β,γ)

Solve the equation (αα'-ββ')x + γ = 0.

Example

julia> A = [3. 4.; 5. 6.]
 2×2 Array{Float64,2}:
  3.0  4.0
  5.0  6.0
@@ -81,4 +81,4 @@
 julia> A*X*A' - E*X*E' + C
 2×2 Array{Float64,2}:
  -2.22045e-16  -4.44089e-16
- -1.33227e-15   1.11022e-15
source
MatrixEquations.lyapds! — Function
lyapds!(A, C; adj = false)

Solve the discrete Lyapunov matrix equation

            op(A)Xop(A)' - X + C = 0,

where op(A) = A if adj = false and op(A) = A' if adj = true. A is in a real or complex Schur form and C is a symmetric or hermitian matrix. A must not have two eigenvalues α and β such that αβ = 1. The computed symmetric or hermitian solution X is contained in C.

source
lyapds!(A, E, C; adj = false)

Solve the generalized discrete Lyapunov matrix equation

            op(A)Xop(A)' - op(E)Xop(E)' + C = 0,

where op(A) = A and op(E) = E if adj = false and op(A) = A' and op(E) = E' if adj = true. The pair (A,E) is in a generalized real or complex Schur form and C is a symmetric or hermitian matrix. The pencil A-λE must not have two eigenvalues α and β such that αβ = 1. The computed symmetric or hermitian solution X is contained in C.

source
+ -1.33227e-15 1.11022e-15source
MatrixEquations.lyapds! — Function
lyapds!(A, C; adj = false)

Solve the discrete Lyapunov matrix equation

            op(A)Xop(A)' - X + C = 0,

where op(A) = A if adj = false and op(A) = A' if adj = true. A is in a real or complex Schur form and C is a symmetric or hermitian matrix. A must not have two eigenvalues α and β such that αβ = 1. The computed symmetric or hermitian solution X is contained in C.

source
lyapds!(A, E, C; adj = false)

Solve the generalized discrete Lyapunov matrix equation

            op(A)Xop(A)' - op(E)Xop(E)' + C = 0,

where op(A) = A and op(E) = E if adj = false and op(A) = A' and op(E) = E' if adj = true. The pair (A,E) is in a generalized real or complex Schur form and C is a symmetric or hermitian matrix. The pencil A-λE must not have two eigenvalues α and β such that αβ = 1. The computed symmetric or hermitian solution X is contained in C.

source
diff --git a/dev/makeindex.html b/dev/makeindex.html index 29c46f5..34c2349 100644 --- a/dev/makeindex.html +++ b/dev/makeindex.html @@ -1,2 +1,2 @@ -Index · MatrixEquations.jl
GitHub

Index

+Index · MatrixEquations.jl
GitHub

Index

diff --git a/dev/meoperators.html b/dev/meoperators.html index 0be231c..4241e90 100644 --- a/dev/meoperators.html +++ b/dev/meoperators.html @@ -1,5 +1,5 @@ -Linear Operators Related to Matrix Equation Solvers · MatrixEquations.jl
GitHub

Linear Operators Related to Matrix Equation Solvers

Lyapunov and Lyapunov-like Operators

MatrixEquations.lyapop — Function
L = lyapop(A; disc = false, her = false)

Define, for an n x n matrix A, the continuous Lyapunov operator L:X -> AX+XA' if disc = false or the discrete Lyapunov operator L:X -> AXA'-X if disc = true. If her = false the Lyapunov operator L:X -> Y maps general square matrices X into general square matrices Y, and the associated matrix M = Matrix(L) is $n^2 \times n^2$. If her = true the Lyapunov operator L:X -> Y maps symmetric/Hermitian matrices X into symmetric/Hermitian matrices Y, and the associated matrix M = Matrix(L) is $n(n+1)/2 \times n(n+1)/2$. For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov operators. Linear Algebra and its Applications 312:35–71, 2000.

source
L = lyapop(A, E; disc = false, her = false)

Define, for a pair (A,E) of n x n matrices, the continuous Lyapunov operator L:X -> AXE'+EXA' if disc = false or the discrete Lyapunov operator L:X -> AXA'-EXE' if disc = true. If her = false the Lyapunov operator L:X -> Y maps general square matrices X into general square matrices Y, and the associated matrix M = Matrix(L) is $n^2 \times n^2$. If her = true the Lyapunov operator L:X -> Y maps symmetric/Hermitian matrices X into symmetric/Hermitian matrices Y, and the associated M = Matrix(L) is a $n(n+1)/2 \times n(n+1)/2$. For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov operators. Linear Algebra and its Applications 312:35–71, 2000.

source
MatrixEquations.invlyapop — Function
LINV = invlyapop(A; disc = false, her = false)

Define LINV, the inverse of the continuous Lyapunov operator L:X -> AX+XA' for disc = false or the inverse of the discrete Lyapunov operator L:X -> AXA'-X for disc = true, where A is an n x n matrix. If her = false the inverse Lyapunov operator LINV:Y -> X maps general square matrices Y into general square matrices X, and the associated matrix M = Matrix(LINV) is $n^2 \times n^2$. If her = true the inverse Lyapunov operator LINV:Y -> X maps symmetric/Hermitian matrices Y into symmetric/Hermitian matrices X, and the associated matrix M = Matrix(LINV) is $n(n+1)/2 \times n(n+1)/2$. For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov operators. Linear Algebra and its Applications 312:35–71, 2000.

source
LINV = invlyapop(A, E; disc = false, her = false)

Define LINV, the inverse of the continuous Lyapunov operator L:X -> AXE'+EXA' for disc = false or the inverse of the discrete Lyapunov operator L:X -> AXA'-EXE' for disc = true, where (A,E) is a pair of n x n matrices. If her = false the inverse Lyapunov operator LINV:Y -> X maps general square matrices Y into general square matrices X, and the associated matrix M = Matrix(LINV) is $n^2 \times n^2$. If her = true the inverse Lyapunov operator LINV:Y -> X maps symmetric/Hermitian matrices Y into symmetric/Hermitian matrices X, and the associated matrix M = Matrix(LINV) is $n(n+1)/2 \times n(n+1)/2$. For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov operators. Linear Algebra and its Applications 312:35–71, 2000.

source
MatrixEquations.lyaplikeop — Function
L = lyaplikeop(A; isig = 1, adj = false, htype = false)

For a matrix A, define for adj = false the continuous T-Lyapunov operator L:X -> A*X+isig*transpose(X)*transpose(A) if htype = false or the continuous H-Lyapunov operator L:X -> A*X+isig*X'*A' if htype = true, or define for adj = true the continuous T-Lyapunov operator L:X -> A*transpose(X)+X*transpose(A) if htype = false, or the continuous H-Lyapunov operator L:X -> A*X'+isig*X*A' if htype = true.

source
MatrixEquations.tulyaplikeop — Function
L = tulyaplikeop(U; adj = false)

Define, for an upper triangular matrix U, the continuous T-Lyapunov operator L:X -> transpose(U)*X+transpose(X)*U, if adj = false, or L:X -> U*transpose(X)+X*transpose(U) if adj = true.

source
MatrixEquations.hulyaplikeop — Function
L = hulyaplikeop(U; adj = false)

Define, for an upper triangular matrix U, the continuous H-Lyapunov operator L:X -> U'*X+X'*U, if adj = false, L:X -> U*X'+X*U' if adj = true.

source

Sylvester and Sylvester-like Operators

MatrixEquations.sylvop — Function
M = sylvop(A, B; disc = false)

Define the continuous Sylvester operator M: X -> AX+XB if disc = false or the discrete Sylvester operator M: X -> AXB+X if disc = true, where A and B are square matrices.

source
M = sylvop(A, B, C, D)

Define the generalized Sylvester operator M: X -> AXB+CXD, where (A,C) and (B,D) are pairs of square matrices.

source
MatrixEquations.invsylvop — Function
MINV = invsylvop(A, B; disc = false)

Define MINV, the inverse of the continuous Sylvester operator M: X -> AX+XB if disc = false or of the discrete Sylvester operator M: X -> AXB+X if disc = true, where A and B are square matrices.

source
MINV = invsylvop(A, B, C, D)

Define MINV, the inverse of the generalized Sylvester operator M: X -> AXB+CXD, where (A,C) and (B,D) are pairs of square matrices.

source
MatrixEquations.gsylvop — Function
M = gsylvop(A, B, C, D; mx, nx, htype = false)

Define the generalized T-Sylvester operator M: X -> ∑ A_i*X*B_i + ∑ C_j'*transpose(X)*D_j, if htype = false or the generalized H-Sylvester operator M: X -> ∑ A_i*X*B_i + ∑ C_j'*X'*D_j, if htype = true. A_i and C_j are matrices having the same row dimension and B_i and D_j are matrices having the same column dimension. A_i and B_i are contained in the k-vectors of matrices A and B, respectively, and C_j and D_j are contained in the l-vectors of matrices C and D, respectively. Any of the component matrices can be given as an UniformScaling. The keyword parameters mx and nx can be used to specify the row and column dimensions of X, if they cannot be inferred from the data.

source

Sylvester System Operators

MatrixEquations.sylvsysop — Function
M = sylvsysop(A, B, C, D)

Define the operator M: (X,Y) -> (AX+YB, CX+YD), where (A,C) and (B,D) are pairs of square matrices.

source
MatrixEquations.invsylvsysop — Function
MINV = invsylvsysop(A, B, C, D)

Define MINV, the inverse of the linear operator M: (X,Y) -> (AX+YB, CX+YD ), where (A,C) and (B,D) a pairs of square matrices.

source

Elementary Matrix Operators

MatrixEquations.trmatop — Type
M = trmatop(n, m)
-M = trmatop(A)

Define the transposition operator M: X -> X' for all n x m matrices or for all matrices of size of A. The corresponding commutation matrix (see here) can be generated as Matrix(M).

source
MatrixEquations.eliminationop — Type
M = eliminationop(n)
-M = eliminationop(A)

Define the elimination operator of all n×n matrices to select their upper triangular parts or of all square matrices of size of A. See here for the definition of an elimination matrix for the selection of lower triangular parts.

source
MatrixEquations.duplicationop — Type
M = duplicationop(n)
-M = duplicationop(A)

Define the duplication operator of all n×n matrices to reconstruct a hermitian matrix from its upper triangular elements or of all square matrices of size of A. See here for the definition of a duplication matrix from the lower triangular parts.

source
+Linear Operators Related to Matrix Equation Solvers · MatrixEquations.jl
GitHub

Linear Operators Related to Matrix Equation Solvers

Lyapunov and Lyapunov-like Operators

MatrixEquations.lyapop — Function
L = lyapop(A; disc = false, her = false)

Define, for an n x n matrix A, the continuous Lyapunov operator L:X -> AX+XA' if disc = false or the discrete Lyapunov operator L:X -> AXA'-X if disc = true. If her = false the Lyapunov operator L:X -> Y maps general square matrices X into general square matrices Y, and the associated matrix M = Matrix(L) is $n^2 \times n^2$. If her = true the Lyapunov operator L:X -> Y maps symmetric/Hermitian matrices X into symmetric/Hermitian matrices Y, and the associated matrix M = Matrix(L) is $n(n+1)/2 \times n(n+1)/2$. For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov operators. Linear Algebra and its Applications 312:35–71, 2000.

source
L = lyapop(A, E; disc = false, her = false)

Define, for a pair (A,E) of n x n matrices, the continuous Lyapunov operator L:X -> AXE'+EXA' if disc = false or the discrete Lyapunov operator L:X -> AXA'-EXE' if disc = true. If her = false the Lyapunov operator L:X -> Y maps general square matrices X into general square matrices Y, and the associated matrix M = Matrix(L) is $n^2 \times n^2$. If her = true the Lyapunov operator L:X -> Y maps symmetric/Hermitian matrices X into symmetric/Hermitian matrices Y, and the associated M = Matrix(L) is a $n(n+1)/2 \times n(n+1)/2$. For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov operators. Linear Algebra and its Applications 312:35–71, 2000.

source
MatrixEquations.invlyapop — Function
LINV = invlyapop(A; disc = false, her = false)

Define LINV, the inverse of the continuous Lyapunov operator L:X -> AX+XA' for disc = false or the inverse of the discrete Lyapunov operator L:X -> AXA'-X for disc = true, where A is an n x n matrix. If her = false the inverse Lyapunov operator LINV:Y -> X maps general square matrices Y into general square matrices X, and the associated matrix M = Matrix(LINV) is $n^2 \times n^2$. If her = true the inverse Lyapunov operator LINV:Y -> X maps symmetric/Hermitian matrices Y into symmetric/Hermitian matrices X, and the associated matrix M = Matrix(LINV) is $n(n+1)/2 \times n(n+1)/2$. For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov operators. Linear Algebra and its Applications 312:35–71, 2000.

source
LINV = invlyapop(A, E; disc = false, her = false)

Define LINV, the inverse of the continuous Lyapunov operator L:X -> AXE'+EXA' for disc = false or the inverse of the discrete Lyapunov operator L:X -> AXA'-EXE' for disc = true, where (A,E) is a pair of n x n matrices. If her = false the inverse Lyapunov operator LINV:Y -> X maps general square matrices Y into general square matrices X, and the associated matrix M = Matrix(LINV) is $n^2 \times n^2$. If her = true the inverse Lyapunov operator LINV:Y -> X maps symmetric/Hermitian matrices Y into symmetric/Hermitian matrices X, and the associated matrix M = Matrix(LINV) is $n(n+1)/2 \times n(n+1)/2$. For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov operators. Linear Algebra and its Applications 312:35–71, 2000.

source
MatrixEquations.lyaplikeop — Function
L = lyaplikeop(A; isig = 1, adj = false, htype = false)

For a matrix A, define for adj = false the continuous T-Lyapunov operator L:X -> A*X+isig*transpose(X)*transpose(A) if htype = false or the continuous H-Lyapunov operator L:X -> A*X+isig*X'*A' if htype = true, or define for adj = true the continuous T-Lyapunov operator L:X -> A*transpose(X)+X*transpose(A) if htype = false, or the continuous H-Lyapunov operator L:X -> A*X'+isig*X*A' if htype = true.

source
MatrixEquations.tulyaplikeop — Function
L = tulyaplikeop(U; adj = false)

Define, for an upper triangular matrix U, the continuous T-Lyapunov operator L:X -> transpose(U)*X+transpose(X)*U, if adj = false, or L:X -> U*transpose(X)+X*transpose(U) if adj = true.

source
MatrixEquations.hulyaplikeop — Function
L = hulyaplikeop(U; adj = false)

Define, for an upper triangular matrix U, the continuous H-Lyapunov operator L:X -> U'*X+X'*U, if adj = false, L:X -> U*X'+X*U' if adj = true.

source

Sylvester and Sylvester-like Operators

MatrixEquations.sylvop — Function
M = sylvop(A, B; disc = false)

Define the continuous Sylvester operator M: X -> AX+XB if disc = false or the discrete Sylvester operator M: X -> AXB+X if disc = true, where A and B are square matrices.

source
M = sylvop(A, B, C, D)

Define the generalized Sylvester operator M: X -> AXB+CXD, where (A,C) and (B,D) are pairs of square matrices.

source
MatrixEquations.invsylvop — Function
MINV = invsylvop(A, B; disc = false)

Define MINV, the inverse of the continuous Sylvester operator M: X -> AX+XB if disc = false or of the discrete Sylvester operator M: X -> AXB+X if disc = true, where A and B are square matrices.

source
MINV = invsylvop(A, B, C, D)

Define MINV, the inverse of the generalized Sylvester operator M: X -> AXB+CXD, where (A,C) and (B,D) are pairs of square matrices.

source
MatrixEquations.gsylvop — Function
M = gsylvop(A, B, C, D; mx, nx, htype = false)

Define the generalized T-Sylvester operator M: X -> ∑ A_i*X*B_i + ∑ C_j'*transpose(X)*D_j, if htype = false or the generalized H-Sylvester operator M: X -> ∑ A_i*X*B_i + ∑ C_j'*X'*D_j, if htype = true. A_i and C_j are matrices having the same row dimension and B_i and D_j are matrices having the same column dimension. A_i and B_i are contained in the k-vectors of matrices A and B, respectively, and C_j and D_j are contained in the l-vectors of matrices C and D, respectively. Any of the component matrices can be given as an UniformScaling. The keyword parameters mx and nx can be used to specify the row and column dimensions of X, if they cannot be inferred from the data.

source

Sylvester System Operators

MatrixEquations.sylvsysop — Function
M = sylvsysop(A, B, C, D)

Define the operator M: (X,Y) -> (AX+YB, CX+YD), where (A,C) and (B,D) are pairs of square matrices.

source
MatrixEquations.invsylvsysop — Function
MINV = invsylvsysop(A, B, C, D)

Define MINV, the inverse of the linear operator M: (X,Y) -> (AX+YB, CX+YD ), where (A,C) and (B,D) a pairs of square matrices.

source

Elementary Matrix Operators

MatrixEquations.trmatop — Type
M = trmatop(n, m)
+M = trmatop(A)

Define the transposition operator M: X -> X' for all n x m matrices or for all matrices of size of A. The corresponding commutation matrix (see here) can be generated as Matrix(M).

source
MatrixEquations.eliminationop — Type
M = eliminationop(n)
+M = eliminationop(A)

Define the elimination operator of all n×n matrices to select their upper triangular parts or of all square matrices of size of A. See here for the definition of an elimination matrix for the selection of lower triangular parts.

source
MatrixEquations.duplicationop — Type
M = duplicationop(n)
+M = duplicationop(A)

Define the duplication operator of all n×n matrices to reconstruct a hermitian matrix from its upper triangular elements or of all square matrices of size of A. See here for the definition of a duplication matrix from the lower triangular parts.

source
diff --git a/dev/meutil.html b/dev/meutil.html index 6e62a30..1b85679 100644 --- a/dev/meutil.html +++ b/dev/meutil.html @@ -1,2 +1,2 @@ -Matrix Equations Utilities · MatrixEquations.jl
GitHub

Matrix Equations Utilities

MatrixEquations.utqu — Function
X = utqu(Q,U)

Compute efficiently the symmetric/hermitian product X = U'QU, where Q is a symmetric/hermitian matrix.

source
MatrixEquations.utqu! — Function
utqu!(Q,U) -> Q

Compute efficiently the symmetric/hermitian product U'QU -> Q, where Q is a symmetric/hermitian matrix and U is a square matrix. The resulting product overwrites Q.

source
MatrixEquations.qrupdate! — Function
qrupdate!(R, Y) -> R

Update the upper triangular factor R by the upper triangular factor of the QR factorization of [ R; Y' ], where Y is a low-rank matrix Y (typically with one or two columns). The computation of R only uses O(n^2) operations (n is the size of R). The input matrix R is updated in place and the matrix Y is destroyed during the computation.

source
MatrixEquations.rqupdate! — Function
rqupdate!(R, Y) -> R

Update the upper triangular factor R by the upper triangular factor of the RQ factorization of [ Y R], where Y is a low-rank matrix Y (typically with one or two columns). The computation of R only uses O(n^2) operations (n is the size of R). The input matrix R is updated in place and the matrix Y is destroyed during the computation.

source
MatrixEquations.isschur — Function
 isschur(A::AbstractMatrix) -> Bool

Test whether A is a square matrix in a real or complex Schur form. In the real case, it is only tested whether A is a quasi upper triangular matrix, which may have 2x2 diagonal blocks, which however must not correspond to complex conjugate eigenvalues. In the complex case, it is tested if A is upper triangular.

source
 isschur(A::AbstractMatrix, B::AbstractMatrix) -> Bool

Test whether (A,B) is a pair of square matrices in a real or complex generalized Schur form. In the real case, it is only tested whether B is upper triangular and A is a quasi upper triangular matrix, which may have 2x2 diagonal blocks, which however must not correspond to complex conjugate eigenvalues. In the complex case, it is tested if A and B are both upper triangular.

source
MatrixEquations.triu2vec — Function
x = triu2vec(Q; rowwise = false, her = false)

Reshape the upper triangular part of the nxn array Q as a one-dimensional column vector x with n(n+1)/2 elements. Q is assumed symmetric/hermitian if her = true. The elements of x correspond to stacking the elements of successive columns of the upper triangular part of Q, if rowwise = false, or stacking the elements of successive rows of the upper triangular part of Q, if rowwise = true.

source
MatrixEquations.vec2triu — Function
Q = vec2triu(x; rowwise = false, her = false)

Build from a one-dimensional column vector x with n(n+1)/2 elements an nxn upper triangular matrix Q if her = false or an nxn symmetric/hermitian array Q if her = true. The elements of x correspond to stacking the elements of successive columns of the upper triangular part of Q, if rowwise = false, or stacking the elements of successive rows of the upper triangular part of Q, if rowwise = true.

source
+Matrix Equations Utilities · MatrixEquations.jl
GitHub

Matrix Equations Utilities

MatrixEquations.utqu — Function
X = utqu(Q,U)

Compute efficiently the symmetric/hermitian product X = U'QU, where Q is a symmetric/hermitian matrix.

source
MatrixEquations.utqu! — Function
utqu!(Q,U) -> Q

Compute efficiently the symmetric/hermitian product U'QU -> Q, where Q is a symmetric/hermitian matrix and U is a square matrix. The resulting product overwrites Q.

source
MatrixEquations.qrupdate! — Function
qrupdate!(R, Y) -> R

Update the upper triangular factor R by the upper triangular factor of the QR factorization of [ R; Y' ], where Y is a low-rank matrix Y (typically with one or two columns). The computation of R only uses O(n^2) operations (n is the size of R). The input matrix R is updated in place and the matrix Y is destroyed during the computation.

source
MatrixEquations.rqupdate! — Function
rqupdate!(R, Y) -> R

Update the upper triangular factor R by the upper triangular factor of the RQ factorization of [ Y R], where Y is a low-rank matrix Y (typically with one or two columns). The computation of R only uses O(n^2) operations (n is the size of R). The input matrix R is updated in place and the matrix Y is destroyed during the computation.

source
MatrixEquations.isschur — Function
 isschur(A::AbstractMatrix) -> Bool

Test whether A is a square matrix in a real or complex Schur form. In the real case, it is only tested whether A is a quasi upper triangular matrix, which may have 2x2 diagonal blocks, which however must not correspond to complex conjugate eigenvalues. In the complex case, it is tested if A is upper triangular.

source
 isschur(A::AbstractMatrix, B::AbstractMatrix) -> Bool

Test whether (A,B) is a pair of square matrices in a real or complex generalized Schur form. In the real case, it is only tested whether B is upper triangular and A is a quasi upper triangular matrix, which may have 2x2 diagonal blocks, which however must not correspond to complex conjugate eigenvalues. In the complex case, it is tested if A and B are both upper triangular.

source
MatrixEquations.triu2vec — Function
x = triu2vec(Q; rowwise = false, her = false)

Reshape the upper triangular part of the nxn array Q as a one-dimensional column vector x with n(n+1)/2 elements. Q is assumed symmetric/hermitian if her = true. The elements of x correspond to stacking the elements of successive columns of the upper triangular part of Q, if rowwise = false, or stacking the elements of successive rows of the upper triangular part of Q, if rowwise = true.

source
MatrixEquations.vec2triu — Function
Q = vec2triu(x; rowwise = false, her = false)

Build from a one-dimensional column vector x with n(n+1)/2 elements an nxn upper triangular matrix Q if her = false or an nxn symmetric/hermitian array Q if her = true. The elements of x correspond to stacking the elements of successive columns of the upper triangular part of Q, if rowwise = false, or stacking the elements of successive rows of the upper triangular part of Q, if rowwise = true.

source
diff --git a/dev/plyapunov.html b/dev/plyapunov.html index 3153a27..3a67244 100644 --- a/dev/plyapunov.html +++ b/dev/plyapunov.html @@ -19,7 +19,7 @@ julia> A*U*U'+U*U'*A'+B*B' 2×2 Array{Float64,2}: 0.0 8.88178e-16 - 8.88178e-16 3.55271e-15source
U = plyapc(A, E, B)

Compute U, the upper triangular factor of the solution X = UU' of the generalized continuous Lyapunov equation

  AXE' + EXA' + BB' = 0,

where A and E are square real or complex matrices and B is a matrix with the same number of rows as A. The pencil A - λE must have only eigenvalues with negative real parts.

U = plyapc(A', E', B')

Compute U, the upper triangular factor of the solution X = U'U of the generalized continuous Lyapunov equation

  A'XE + E'XA + B'B = 0,

where A and E are square real or complex matrices and B is a matrix with the same number of columns as A. The pencil A - λE must have only eigenvalues with negative real parts.

Example

julia> using LinearAlgebra
+ 8.88178e-16  3.55271e-15
source
U = plyapc(A, E, B)

Compute U, the upper triangular factor of the solution X = UU' of the generalized continuous Lyapunov equation

  AXE' + EXA' + BB' = 0,

where A and E are square real or complex matrices and B is a matrix with the same number of rows as A. The pencil A - λE must have only eigenvalues with negative real parts.

U = plyapc(A', E', B')

Compute U, the upper triangular factor of the solution X = U'U of the generalized continuous Lyapunov equation

  A'XE + E'XA + B'B = 0,

where A and E are square real or complex matrices and B is a matrix with the same number of columns as A. The pencil A - λE must have only eigenvalues with negative real parts.

Example

julia> using LinearAlgebra
 
 julia> A = [-2. 1.;-1. -2.]
 2×2 Array{Float64,2}:
@@ -44,7 +44,7 @@
 julia> A*U*U'*E'+E*U*U'*A'+B*B'
 2×2 Array{Float64,2}:
   0.0          -8.88178e-16
- -1.33227e-15  -2.66454e-15
source
MatrixEquations.plyapcs! — Function
plyapcs!(A,R;adj = false)

Solve the positive continuous Lyapunov matrix equation

            op(A)X + Xop(A)' + op(R)*op(R)' = 0

for X = op(U)*op(U)', where op(K) = K if adj = false and op(K) = K' if adj = true. A is a square real matrix in a real Schur form or a square complex matrix in a complex Schur form and R is an upper triangular matrix. A must have only eigenvalues with negative real parts. R contains on output the solution U.

source
plyapcs!(A,E,R;adj = false)

Solve the generalized positive continuous Lyapunov matrix equation

            op(A)Xop(E)' + op(E)*Xop(A)' + op(R)*op(R)' = 0

for X = op(U)*op(U)', where op(K) = K if adj = false and op(K) = K' if adj = true. The pair (A,E) is in a generalized real/complex Schur form and R is an upper triangular matrix. The pencil A-λE must have only eigenvalues with negative real parts. R contains on output the solution U.

source

Discrete-time Lyapunov (Stein) Matrix Equations

MatrixEquations.plyapd — Function
U = plyapd(A, B)

Compute U, the upper triangular factor of the solution X = UU' of the discrete Lyapunov equation

  AXA' - X + BB' = 0,

where A is a square real or complex matrix and B is a matrix with the same number of rows as A. A must have only eigenvalues with moduli less than one.

U = plyapd(A', B')

Compute U, the upper triangular factor of the solution X = U'U of the discrete Lyapunov equation

  A'XA - X + B'B = 0,

where A is a square real or complex matrix and B is a matrix with the same number of columns as A. A must have only eigenvalues with moduli less than one.

Example

julia> using LinearAlgebra
+ -1.33227e-15  -2.66454e-15
source
MatrixEquations.plyapcs! — Function
plyapcs!(A,R;adj = false)

Solve the positive continuous Lyapunov matrix equation

            op(A)X + Xop(A)' + op(R)*op(R)' = 0

for X = op(U)*op(U)', where op(K) = K if adj = false and op(K) = K' if adj = true. A is a square real matrix in a real Schur form or a square complex matrix in a complex Schur form and R is an upper triangular matrix. A must have only eigenvalues with negative real parts. R contains on output the solution U.

source
plyapcs!(A,E,R;adj = false)

Solve the generalized positive continuous Lyapunov matrix equation

            op(A)Xop(E)' + op(E)*Xop(A)' + op(R)*op(R)' = 0

for X = op(U)*op(U)', where op(K) = K if adj = false and op(K) = K' if adj = true. The pair (A,E) is in a generalized real/complex Schur form and R is an upper triangular matrix. The pencil A-λE must have only eigenvalues with negative real parts. R contains on output the solution U.

source

Discrete-time Lyapunov (Stein) Matrix Equations

MatrixEquations.plyapd — Function
U = plyapd(A, B)

Compute U, the upper triangular factor of the solution X = UU' of the discrete Lyapunov equation

  AXA' - X + BB' = 0,

where A is a square real or complex matrix and B is a matrix with the same number of rows as A. A must have only eigenvalues with moduli less than one.

U = plyapd(A', B')

Compute U, the upper triangular factor of the solution X = U'U of the discrete Lyapunov equation

  A'XA - X + B'B = 0,

where A is a square real or complex matrix and B is a matrix with the same number of columns as A. A must have only eigenvalues with moduli less than one.

Example

julia> using LinearAlgebra
 
 julia> A = [-0.5 .1;-0.1 -0.5]
 2×2 Array{Float64,2}:
@@ -64,7 +64,7 @@
 julia> A*U*U'*A'-U*U'+B*B'
 2×2 Array{Float64,2}:
  -4.44089e-16  4.44089e-16
-  4.44089e-16  1.77636e-15
source
U = plyapd(A, E, B)

Compute U, the upper triangular factor of the solution X = UU' of the generalized discrete Lyapunov equation

  AXA' - EXE' + BB' = 0,

where A and E are square real or complex matrices and B is a matrix with the same number of rows as A. The pencil A - λE must have only eigenvalues with moduli less than one.

U = plyapd(A', E', B')

Compute U, the upper triangular factor of the solution X = U'U of the generalized discrete Lyapunov equation

  A'XA - E'XE + B'B = 0,

where A and E are square real or complex matrices and B is a matrix with the same number of columns as A. The pencil A - λE must have only eigenvalues with moduli less than one.

Example

julia> using LinearAlgebra
+  4.44089e-16  1.77636e-15
source
U = plyapd(A, E, B)

Compute U, the upper triangular factor of the solution X = UU' of the generalized discrete Lyapunov equation

  AXA' - EXE' + BB' = 0,

where A and E are square real or complex matrices and B is a matrix with the same number of rows as A. The pencil A - λE must have only eigenvalues with moduli less than one.

U = plyapd(A', E', B')

Compute U, the upper triangular factor of the solution X = U'U of the generalized discrete Lyapunov equation

  A'XA - E'XE + B'B = 0,

where A and E are square real or complex matrices and B is a matrix with the same number of columns as A. The pencil A - λE must have only eigenvalues with moduli less than one.

Example

julia> using LinearAlgebra
 
 julia> A = [-0.5 .1;-0.1 -0.5]
 2×2 Array{Float64,2}:
@@ -89,4 +89,4 @@
 julia> A*U*U'*A'-E*U*U'*E'+B*B'
 2×2 Array{Float64,2}:
  1.77636e-15  2.22045e-15
- 2.22045e-15  2.66454e-15
source
MatrixEquations.plyapds! — Function
plyapds!(A, R; adj = false)

Solve the positive discrete Lyapunov matrix equation

            op(A)Xop(A)' - X + op(R)*op(R)' = 0

for X = op(U)*op(U)', where op(K) = K if adj = false and op(K) = K' if adj = true. A is a square real matrix in a real Schur form or a square complex matrix in a complex Schur form and R is an upper triangular matrix. A must have only eigenvalues with moduli less than one. R contains on output the upper triangular solution U.

source
plyapds!(A,E,R;adj = false)

Solve the generalized positive discrete Lyapunov matrix equation

            op(A)Xop(A)' - op(E)Xop(E)' + op(R)*op(R)' = 0

for X = op(U)*op(U)', where op(K) = K if adj = false and op(K) = K' if adj = true. The pair (A,E) of square real or complex matrices is in a generalized Schur form and R is an upper triangular matrix. A-λE must have only eigenvalues with moduli less than one. R contains on output the upper triangular solution U.

source

Schur Form Based Solvers

MatrixEquations.plyaps — Function
U = plyaps(A, B; disc = false)

Compute U, the upper triangular factor of the solution X = UU' of the continuous Lyapunov equation

  AX + XA' + BB' = 0,

where A is a square real or complex matrix in a real or complex Schur form, respectively, and B is a matrix with the same number of rows as A. A must have only eigenvalues with negative real parts. Only the upper Hessenberg part of A is referenced.

U = plyaps(A', B', disc = false)

Compute U, the upper triangular factor of the solution X = U'U of the continuous Lyapunov equation

  A'X + XA + B'B = 0,

where A is a square real or complex matrix in a real or complex Schur form, respectively, and B is a matrix with the same number of columns as A. A must have only eigenvalues with negative real parts. Only the upper Hessenberg part of A is referenced.

U = plyaps(A, B, disc = true)

Compute U, the upper triangular factor of the solution X = UU' of the discrete Lyapunov equation

  AXA' - X + BB' = 0,

where A is a square real or complex matrix in a real or complex Schur form, respectively, and B is a matrix with the same number of rows as A. A must have only eigenvalues with moduli less than one. Only the upper Hessenberg part of A is referenced.

U = plyaps(A', B', disc = true)

Compute U, the upper triangular factor of the solution X = U'U of the discrete Lyapunov equation

  A'XA - X + B'B = 0,

where A is a square real or complex matrix in a real or complex Schur form, respectively, and B is a matrix with the same number of columns as A. A must have only eigenvalues with moduli less than one. Only the upper Hessenberg part of A is referenced.

source
U = plyaps(A, E, B; disc = false)

Compute U, the upper triangular factor of the solution X = UU' of the generalized continuous Lyapunov equation

  AXE' + EXA' + BB' = 0,

where A and E are square real or complex matrices with the pair (A,E) in a generalied real or complex Schur form, respectively, and B is a matrix with the same number of rows as A. The pencil A - λE must have only eigenvalues with negative real parts.

U = plyaps(A', E', B', disc = false)

Compute U, the upper triangular factor of the solution X = U'U of the generalized continuous Lyapunov equation

  A'XE + E'XA + B'B = 0,

where A and E are square real or complex matrices with the pair (A,E) in a generalied real or complex Schur form, respectively, and B is a matrix with the same number of columns as A. The pencil A - λE must have only eigenvalues with negative real parts.

U = plyaps(A, E, B, disc = true)

Compute U, the upper triangular factor of the solution X = UU' of the generalized discrete Lyapunov equation

  AXA' - EXE' + BB' = 0,

where A and E are square real or complex matrices with the pair (A,E) in a generalied real or complex Schur form, respectively, and B is a matrix with the same number of rows as A. The pencil A - λE must have only eigenvalues with moduli less than one.

U = plyaps(A', E', B', disc = true)

Compute U, the upper triangular factor of the solution X = U'U of the generalized discrete Lyapunov equation

  A'XA - E'XE + B'B = 0,

where A and E are square real or complex matrices with the pair (A,E) in a generalied real or complex Schur form, respectively, and B is a matrix with the same number of columns as A. The pencil A - λE must have only eigenvalues with moduli less than one.

source
+ 2.22045e-15 2.66454e-15source
MatrixEquations.plyapds! — Function
plyapds!(A, R; adj = false)

Solve the positive discrete Lyapunov matrix equation

            op(A)Xop(A)' - X + op(R)*op(R)' = 0

for X = op(U)*op(U)', where op(K) = K if adj = false and op(K) = K' if adj = true. A is a square real matrix in a real Schur form or a square complex matrix in a complex Schur form and R is an upper triangular matrix. A must have only eigenvalues with moduli less than one. R contains on output the upper triangular solution U.

source
plyapds!(A,E,R;adj = false)

Solve the generalized positive discrete Lyapunov matrix equation

            op(A)Xop(A)' - op(E)Xop(E)' + op(R)*op(R)' = 0

for X = op(U)*op(U)', where op(K) = K if adj = false and op(K) = K' if adj = true. The pair (A,E) of square real or complex matrices is in a generalized Schur form and R is an upper triangular matrix. A-λE must have only eigenvalues with moduli less than one. R contains on output the upper triangular solution U.

source

Schur Form Based Solvers

MatrixEquations.plyaps — Function
U = plyaps(A, B; disc = false)

Compute U, the upper triangular factor of the solution X = UU' of the continuous Lyapunov equation

  AX + XA' + BB' = 0,

where A is a square real or complex matrix in a real or complex Schur form, respectively, and B is a matrix with the same number of rows as A. A must have only eigenvalues with negative real parts. Only the upper Hessenberg part of A is referenced.

U = plyaps(A', B', disc = false)

Compute U, the upper triangular factor of the solution X = U'U of the continuous Lyapunov equation

  A'X + XA + B'B = 0,

where A is a square real or complex matrix in a real or complex Schur form, respectively, and B is a matrix with the same number of columns as A. A must have only eigenvalues with negative real parts. Only the upper Hessenberg part of A is referenced.

U = plyaps(A, B, disc = true)

Compute U, the upper triangular factor of the solution X = UU' of the discrete Lyapunov equation

  AXA' - X + BB' = 0,

where A is a square real or complex matrix in a real or complex Schur form, respectively, and B is a matrix with the same number of rows as A. A must have only eigenvalues with moduli less than one. Only the upper Hessenberg part of A is referenced.

U = plyaps(A', B', disc = true)

Compute U, the upper triangular factor of the solution X = U'U of the discrete Lyapunov equation

  A'XA - X + B'B = 0,

where A is a square real or complex matrix in a real or complex Schur form, respectively, and B is a matrix with the same number of columns as A. A must have only eigenvalues with moduli less than one. Only the upper Hessenberg part of A is referenced.

source
U = plyaps(A, E, B; disc = false)

Compute U, the upper triangular factor of the solution X = UU' of the generalized continuous Lyapunov equation

  AXE' + EXA' + BB' = 0,

where A and E are square real or complex matrices with the pair (A,E) in a generalied real or complex Schur form, respectively, and B is a matrix with the same number of rows as A. The pencil A - λE must have only eigenvalues with negative real parts.

U = plyaps(A', E', B', disc = false)

Compute U, the upper triangular factor of the solution X = U'U of the generalized continuous Lyapunov equation

  A'XE + E'XA + B'B = 0,

where A and E are square real or complex matrices with the pair (A,E) in a generalied real or complex Schur form, respectively, and B is a matrix with the same number of columns as A. The pencil A - λE must have only eigenvalues with negative real parts.

U = plyaps(A, E, B, disc = true)

Compute U, the upper triangular factor of the solution X = UU' of the generalized discrete Lyapunov equation

  AXA' - EXE' + BB' = 0,

where A and E are square real or complex matrices with the pair (A,E) in a generalied real or complex Schur form, respectively, and B is a matrix with the same number of rows as A. The pencil A - λE must have only eigenvalues with moduli less than one.

U = plyaps(A', E', B', disc = true)

Compute U, the upper triangular factor of the solution X = U'U of the generalized discrete Lyapunov equation

  A'XA - E'XE + B'B = 0,

where A and E are square real or complex matrices with the pair (A,E) in a generalied real or complex Schur form, respectively, and B is a matrix with the same number of columns as A. The pencil A - λE must have only eigenvalues with moduli less than one.

source
diff --git a/dev/riccati.html b/dev/riccati.html index b690473..3a0c2e2 100644 --- a/dev/riccati.html +++ b/dev/riccati.html @@ -37,7 +37,7 @@ 3-element Array{Complex{Float64},1}: -4.4115475922960075 - 2.4222082620381076im -4.4115475922960075 + 2.4222082620381076im - -4.337128244724374 + 0.0imsource
arec(A, B, R, Q, S; scaling = 'B', as = false, rtol::Real = nϵ, orth = false) -> (X, EVALS, F, Z)

Computes X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the continuous-time algebraic Riccati equation

 A'X + XA - (XB+S)R^(-1)(B'X+S') + Q = 0,

where R and Q are hermitian/symmetric matrices or uniform scaling operators such that R is nonsingular. Scalar-valued R and Q are interpreted as appropriately sized uniform scaling operators R*I and Q*I. S, if not specified, is set to S = zeros(size(B)). The Schur method of [1] is used.

To enhance the accuracy of computations, a block scaling of matrices R, Q and S is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > norm(B)^2/norm(R). Alternative scaling can be performed using the options scaling = 'S', for a special structure preserving scaling, andscaling = 'G', for a general eigenvalue computation oriented scaling. Experimentally, if orth = true, two scaling procedures can be activated using the options scaling = 'D' andscaling = 'T'. Scaling can be disabled with the choice `scaling = 'N'.

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) eigenvalues of A-BF. F is the stabilizing or anti-stabilizing gain matrix F = R^(-1)(B'X+S'). Z = [ U; V; W ] is a basis for the relevant stable/anti-stable deflating subspace such that X = V/U and F = -W/U. An orthogonal basis Z can be determined, with an increased computational cost, by setting orth = true.

Reference:

[1] Laub, A.J., A Schur Method for Solving Algebraic Riccati equations. IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.

Example

julia> using LinearAlgebra
+  -4.337128244724374 + 0.0im
source
arec(A, B, R, Q, S; scaling = 'B', as = false, rtol::Real = nϵ, orth = false) -> (X, EVALS, F, Z)

Computes X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the continuous-time algebraic Riccati equation

 A'X + XA - (XB+S)R^(-1)(B'X+S') + Q = 0,

where R and Q are hermitian/symmetric matrices or uniform scaling operators such that R is nonsingular. Scalar-valued R and Q are interpreted as appropriately sized uniform scaling operators R*I and Q*I. S, if not specified, is set to S = zeros(size(B)). The Schur method of [1] is used.

To enhance the accuracy of computations, a block scaling of matrices R, Q and S is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > norm(B)^2/norm(R). Alternative scaling can be performed using the options scaling = 'S', for a special structure preserving scaling, andscaling = 'G', for a general eigenvalue computation oriented scaling. Experimentally, if orth = true, two scaling procedures can be activated using the options scaling = 'D' andscaling = 'T'. Scaling can be disabled with the choice `scaling = 'N'.

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) eigenvalues of A-BF. F is the stabilizing or anti-stabilizing gain matrix F = R^(-1)(B'X+S'). Z = [ U; V; W ] is a basis for the relevant stable/anti-stable deflating subspace such that X = V/U and F = -W/U. An orthogonal basis Z can be determined, with an increased computational cost, by setting orth = true.

Reference:

[1] Laub, A.J., A Schur Method for Solving Algebraic Riccati equations. IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.

Example

julia> using LinearAlgebra
 
 julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
 3×3 Array{Float64,2}:
@@ -80,7 +80,7 @@
 3-element Array{Complex{Float64},1}:
   -4.377036283999118 - 2.8107164873731234im
   -4.377036283999118 + 2.8107164873731234im
- -1.8663764577096063 + 0.0im
source
arec(A, B, G, R, Q, S; as = false, rtol::Real = nϵ, orth = false) -> (X, EVALS, F, Z)

Computes X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the continuous-time algebraic Riccati equation

 A'X + XA - XGX - (XB+S)R^(-1)(B'X+S') + Q = 0,

where G, R and Q are hermitian/symmetric matrices or uniform scaling operators such that R is nonsingular. Scalar-valued G, R and Q are interpreted as appropriately sized uniform scaling operators G*I, R*I and Q*I.

To enhance the accuracy of computations, a block oriented scaling of matrices G, Q, R and S is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > max(norm(G), norm(B)^2/norm(R)). Alternative scaling can be performed using the options scaling = 'S', for a special structure preserving scaling, andscaling = 'G', for a general eigenvalue computation oriented scaling. If orth = true, two experimental scaling procedures can be activated using the options scaling = 'D' andscaling = 'T'. Scaling can be disabled with the choice `scaling = 'N'.

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) eigenvalues of A-BF-GX. F is the stabilizing or anti-stabilizing gain matrix F = R^(-1)(B'X+S'). Z = [ U; V; W ] is a basis for the relevant stable/anti-stable deflating subspace such that X = V/U and F = -W/U. An orthogonal basis Z can be determined, with an increased computational cost, by setting orth = true.

Reference: Laub, A.J., A Schur Method for Solving Algebraic Riccati equations. IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.

source
MatrixEquations.ared — Function
ared(A, B, R, Q, S; as = false, rtol::Real = nϵ) -> (X, EVALS, F, Z)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the discrete-time algebraic Riccati equation

A'XA - X - (A'XB+S)(R+B'XB)^(-1)(B'XA+S') + Q = 0,

where R and Q are hermitian/symmetric matrices. Scalar-valued R and Q are interpreted as appropriately sized uniform scaling operators R*I and Q*I. S, if not specified, is set to S = zeros(size(B)).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable) eigenvalues of A-BF. F is the stabilizing gain matrix F = (R+B'XB)^(-1)(B'XA+S'). Z = [ U; V; W ] is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = V/(EU) and F = -W/U.

Reference: W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

Example

julia> using LinearAlgebra
+ -1.8663764577096063 + 0.0im
source
arec(A, B, G, R, Q, S; as = false, rtol::Real = nϵ, orth = false) -> (X, EVALS, F, Z)

Computes X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the continuous-time algebraic Riccati equation

 A'X + XA - XGX - (XB+S)R^(-1)(B'X+S') + Q = 0,

where G, R and Q are hermitian/symmetric matrices or uniform scaling operators such that R is nonsingular. Scalar-valued G, R and Q are interpreted as appropriately sized uniform scaling operators G*I, R*I and Q*I.

To enhance the accuracy of computations, a block oriented scaling of matrices G, Q, R and S is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > max(norm(G), norm(B)^2/norm(R)). Alternative scaling can be performed using the options scaling = 'S', for a special structure preserving scaling, andscaling = 'G', for a general eigenvalue computation oriented scaling. If orth = true, two experimental scaling procedures can be activated using the options scaling = 'D' andscaling = 'T'. Scaling can be disabled with the choice `scaling = 'N'.

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) eigenvalues of A-BF-GX. F is the stabilizing or anti-stabilizing gain matrix F = R^(-1)(B'X+S'). Z = [ U; V; W ] is a basis for the relevant stable/anti-stable deflating subspace such that X = V/U and F = -W/U. An orthogonal basis Z can be determined, with an increased computational cost, by setting orth = true.

Reference: Laub, A.J., A Schur Method for Solving Algebraic Riccati equations. IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.

source
MatrixEquations.ared — Function
ared(A, B, R, Q, S; as = false, rtol::Real = nϵ) -> (X, EVALS, F, Z)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the discrete-time algebraic Riccati equation

A'XA - X - (A'XB+S)(R+B'XB)^(-1)(B'XA+S') + Q = 0,

where R and Q are hermitian/symmetric matrices. Scalar-valued R and Q are interpreted as appropriately sized uniform scaling operators R*I and Q*I. S, if not specified, is set to S = zeros(size(B)).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable) eigenvalues of A-BF. F is the stabilizing gain matrix F = (R+B'XB)^(-1)(B'XA+S'). Z = [ U; V; W ] is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = V/(EU) and F = -W/U.

Reference: W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

Example

julia> using LinearAlgebra
 
 julia> A = [ 0. 1.; 0. 0. ]
 2×2 Array{Float64,2}:
@@ -120,7 +120,7 @@
 julia> eigvals(A-B*F)
 2-element Array{Float64,1}:
  -2.7755575615628914e-16
-  0.5
source

Generalized Riccati Matrix Equations

MatrixEquations.garec — Function
garec(A, E, G, Q = 0; scaling = 'B', as = false, rtol::Real = nϵ) -> (X, EVALS, Z)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the generalized continuous-time algebraic Riccati equation

A'XE + E'XA - E'XGXE + Q = 0,

where G and Q are hermitian/symmetric matrices or uniform scaling operators and E is a nonsingular matrix. Scalar-valued G and Q are interpreted as appropriately sized uniform scaling operators G*I and Q*I.

The generalized Schur method of [1] is used. To enhance the accuracy of computations, a block oriented scaling of matrices G and Q is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > norm(G). Alternative scaling can be performed using the options scaling = 'S', for a special structure preserving scaling, andscaling = 'G', for a general eigenvalue computation oriented scaling. Scaling can be disabled with the choice `scaling = 'N'.

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-GXE,E). Z = [ U; V ] is an orthogonal basis for the stable/anti-stable deflating subspace such that X = V/(EU).

Reference:

[1] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

source
garec(A, E, B, R, Q, S; scaling = 'B', as = false, rtol::Real = nϵ) -> (X, EVALS, F, Z)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the generalized continuous-time algebraic Riccati equation

A'XE + E'XA - (E'XB+S)R^(-1)(B'XE+S') + Q = 0,

where R and Q are hermitian/symmetric matrices such that R is nonsingular, and E is a nonsingular matrix. Scalar-valued R and Q are interpreted as appropriately sized uniform scaling operators R*I and Q*I. S, if not specified, is set to S = zeros(size(B)).

The generalized Schur method of [1] is used. To enhance the accuracy of computations, a block oriented scaling of matrices R, Q and S is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > norm(B)^2/norm(R). Alternative scaling can be performed using the options scaling = 'S', for a special structure preserving scaling, andscaling = 'G', for a general eigenvalue computation oriented scaling. Two experimental scaling procedures can be activated using the options scaling = 'D' andscaling = 'T'. Scaling can be disabled with the choice `scaling = 'N'.

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-BF,E). F is the stabilizing/anti-stabilizing gain matrix F = R^(-1)(B'XE+S'). Z = [ U; V; W ] is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = V/(EU) and F = -W/U.

Reference: W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

Example

julia> using LinearAlgebra
+  0.5
source

Generalized Riccati Matrix Equations

MatrixEquations.garec — Function
garec(A, E, G, Q = 0; scaling = 'B', as = false, rtol::Real = nϵ) -> (X, EVALS, Z)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the generalized continuous-time algebraic Riccati equation

A'XE + E'XA - E'XGXE + Q = 0,

where G and Q are hermitian/symmetric matrices or uniform scaling operators and E is a nonsingular matrix. Scalar-valued G and Q are interpreted as appropriately sized uniform scaling operators G*I and Q*I.

The generalized Schur method of [1] is used. To enhance the accuracy of computations, a block oriented scaling of matrices G and Q is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > norm(G). Alternative scaling can be performed using the options scaling = 'S', for a special structure preserving scaling, andscaling = 'G', for a general eigenvalue computation oriented scaling. Scaling can be disabled with the choice `scaling = 'N'.

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-GXE,E). Z = [ U; V ] is an orthogonal basis for the stable/anti-stable deflating subspace such that X = V/(EU).

Reference:

[1] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

source
garec(A, E, B, R, Q, S; scaling = 'B', as = false, rtol::Real = nϵ) -> (X, EVALS, F, Z)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the generalized continuous-time algebraic Riccati equation

A'XE + E'XA - (E'XB+S)R^(-1)(B'XE+S') + Q = 0,

where R and Q are hermitian/symmetric matrices such that R is nonsingular, and E is a nonsingular matrix. Scalar-valued R and Q are interpreted as appropriately sized uniform scaling operators R*I and Q*I. S, if not specified, is set to S = zeros(size(B)).

The generalized Schur method of [1] is used. To enhance the accuracy of computations, a block oriented scaling of matrices R, Q and S is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > norm(B)^2/norm(R). Alternative scaling can be performed using the options scaling = 'S', for a special structure preserving scaling, andscaling = 'G', for a general eigenvalue computation oriented scaling. Two experimental scaling procedures can be activated using the options scaling = 'D' andscaling = 'T'. Scaling can be disabled with the choice `scaling = 'N'.

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-BF,E). F is the stabilizing/anti-stabilizing gain matrix F = R^(-1)(B'XE+S'). Z = [ U; V; W ] is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = V/(EU) and F = -W/U.

Reference: W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

Example

julia> using LinearAlgebra
 
 julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
 3×3 Array{Float64,2}:
@@ -169,7 +169,7 @@
 3-element Array{Complex{Float64},1}:
  -0.6184265391601462 - 0.29132868445957383im
  -0.6184265391601462 + 0.2913286844595739im
-  -0.216130599644518 + 0.0im
source
garec(A, E, B, G, R, Q, S; scaling = 'B', as = false, rtol::Real = nϵ) -> (X, EVALS, F, Z)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the generalized continuous-time algebraic Riccati equation

A'XE + E'XA - E'XGXE - (E'XB+S)R^(-1)(B'XE+S') + Q = 0,

where G, Q and R are hermitian/symmetric matrices such that R is nonsingular, and E is a nonsingular matrix. Scalar-valued G, R and Q are interpreted as appropriately sized uniform scaling operators G*I, R*I and Q*I.

The generalized Schur method of [1] is used. To enhance the accuracy of computations, a block oriented scaling of matrices G, R, Q and S is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > max(norm(G), norm(B)^2/norm(R)). Alternative scaling can be performed using the options scaling = 'S', for a special structure preserving scaling, andscaling = 'G', for a general eigenvalue computation oriented scaling. Two experimental scaling procedures can be activated using the options scaling = 'D' andscaling = 'T'. Scaling can be disabled with the choice `scaling = 'N'.

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-BF-GXE,E). F is the stabilizing/anti-stabilizing gain matrix F = R^(-1)(B'XE+S'). Z = [ U; V; W ] is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = V/(EU) and F = -W/U.

Reference: W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

source
MatrixEquations.gared — Function
gared(A, E, B, R, Q, S; as = false, rtol::Real = nϵ) -> (X, EVALS, F, Z)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the generalized discrete-time algebraic Riccati equation

A'XA - E'XE - (A'XB+S)(R+B'XB)^(-1)(B'XA+S') + Q = 0,

where R and Q are hermitian/symmetric matrices, and E ist non-singular. Scalar-valued R and Q are interpreted as appropriately sized uniform scaling operators R*I and Q*I. S, if not specified, is set to S = zeros(size(B)).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-BF,E). F is the stabilizing/anti-stabilizing gain matrix F = (R+B'XB)^(-1)(B'XA+S'). Z = [ U; V; W ] is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = V/(EU) and F = -W/U.

Reference: W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

Example

julia> using LinearAlgebra
+  -0.216130599644518 + 0.0im
source
garec(A, E, B, G, R, Q, S; scaling = 'B', as = false, rtol::Real = nϵ) -> (X, EVALS, F, Z)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the generalized continuous-time algebraic Riccati equation

A'XE + E'XA - E'XGXE - (E'XB+S)R^(-1)(B'XE+S') + Q = 0,

where G, Q and R are hermitian/symmetric matrices such that R is nonsingular, and E is a nonsingular matrix. Scalar-valued G, R and Q are interpreted as appropriately sized uniform scaling operators G*I, R*I and Q*I.

The generalized Schur method of [1] is used. To enhance the accuracy of computations, a block oriented scaling of matrices G, R, Q and S is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > max(norm(G), norm(B)^2/norm(R)). Alternative scaling can be performed using the options scaling = 'S', for a special structure preserving scaling, andscaling = 'G', for a general eigenvalue computation oriented scaling. Two experimental scaling procedures can be activated using the options scaling = 'D' andscaling = 'T'. Scaling can be disabled with the choice `scaling = 'N'.

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-BF-GXE,E). F is the stabilizing/anti-stabilizing gain matrix F = R^(-1)(B'XE+S'). Z = [ U; V; W ] is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = V/(EU) and F = -W/U.

Reference: W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

source
MatrixEquations.gared — Function
gared(A, E, B, R, Q, S; as = false, rtol::Real = nϵ) -> (X, EVALS, F, Z)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the generalized discrete-time algebraic Riccati equation

A'XA - E'XE - (A'XB+S)(R+B'XB)^(-1)(B'XA+S') + Q = 0,

where R and Q are hermitian/symmetric matrices, and E ist non-singular. Scalar-valued R and Q are interpreted as appropriately sized uniform scaling operators R*I and Q*I. S, if not specified, is set to S = zeros(size(B)).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-BF,E). F is the stabilizing/anti-stabilizing gain matrix F = (R+B'XB)^(-1)(B'XA+S'). Z = [ U; V; W ] is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = V/(EU) and F = -W/U.

Reference: W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

Example

julia> using LinearAlgebra
 
 julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
 3×3 Array{Float64,2}:
@@ -218,4 +218,4 @@
 3-element Array{Float64,1}:
  -0.5238922629921539
  -0.19053355203423886
- -0.08423561575133902
source
+ -0.08423561575133902source diff --git a/dev/sylvester.html b/dev/sylvester.html index 5bb9ee2..b490dbe 100644 --- a/dev/sylvester.html +++ b/dev/sylvester.html @@ -22,7 +22,7 @@ julia> A*X + X*B - C 2×2 Array{Float64,2}: 2.66454e-15 1.77636e-15 - -3.77476e-15 4.44089e-16source
MatrixEquations.sylvcs! — Function
sylvcs!(A,B,C; adjA = false, adjB = false)

Solve the continuous Sylvester matrix equation

            op(A)X + Xop(B) =  C,

where op(A) = A or op(A) = A' if adjA = false or adjA = true, respectively, and op(B) = B or op(B) = B' if adjB = false or adjB = true, respectively. A and B are square matrices in Schur forms, and A and -B must not have common eigenvalues. C contains on output the solution X.

source
MatrixEquations.sylvd — Function
X = sylvd(A,B,C)

Solve the discrete Sylvester matrix equation

            AXB + X = C

using an extension of the Bartels-Stewart Schur form based approach. A and B are square matrices, and A and -B must not have common reciprocal eigenvalues.

The following particular cases are also adressed:

X = sylvd(α*I,B,C)  or  X = sylvd(α,B,C)

Solve the matrix equation X(αB+I) = C.

X = sylvd(A,β*I,C)   or  X = sylvd(A,β,C)

Solve the matrix equation (βA+I)X = C.

X = sylvd(α*I,β*I,C)  or  X = sylvd(α,β,C)

Solve the matrix equation (αβ+1)X = C.

x = sylvd(α,β,γ)

Solve the equation (αβ+1)x = γ.

Example

julia> A = [3. 4.; 5. 6.]
+ -3.77476e-15  4.44089e-16
source
MatrixEquations.sylvcs! — Function
sylvcs!(A,B,C; adjA = false, adjB = false)

Solve the continuous Sylvester matrix equation

            op(A)X + Xop(B) =  C,

where op(A) = A or op(A) = A' if adjA = false or adjA = true, respectively, and op(B) = B or op(B) = B' if adjB = false or adjB = true, respectively. A and B are square matrices in Schur forms, and A and -B must not have common eigenvalues. C contains on output the solution X.

source
MatrixEquations.sylvd — Function
X = sylvd(A,B,C)

Solve the discrete Sylvester matrix equation

            AXB + X = C

using an extension of the Bartels-Stewart Schur form based approach. A and B are square matrices, and A and -B must not have common reciprocal eigenvalues.

The following particular cases are also adressed:

X = sylvd(α*I,B,C)  or  X = sylvd(α,B,C)

Solve the matrix equation X(αB+I) = C.

X = sylvd(A,β*I,C)   or  X = sylvd(A,β,C)

Solve the matrix equation (βA+I)X = C.

X = sylvd(α*I,β*I,C)  or  X = sylvd(α,β,C)

Solve the matrix equation (αβ+1)X = C.

x = sylvd(α,β,γ)

Solve the equation (αβ+1)x = γ.

Example

julia> A = [3. 4.; 5. 6.]
 2×2 Array{Float64,2}:
  3.0  4.0
  5.0  6.0
@@ -45,7 +45,7 @@
 julia> A*X*B + X - C
 2×2 Array{Float64,2}:
   8.88178e-16   8.88178e-16
- -3.9968e-15   -5.55112e-15
source
MatrixEquations.sylvds! — Function
sylvds!(A,B,C; adjA = false, adjB = false)

Solve the discrete Sylvester matrix equation

            op(A)Xop(B) + X =  C,

where op(A) = A or op(A) = A' if adjA = false or adjA = true, respectively, and op(B) = B or op(B) = B' if adjB = false or adjB = true, respectively. A and B are square matrices in Schur forms, and A and -B must not have common reciprocal eigenvalues. C contains on output the solution X.

source

Generalized Sylvester Matrix Equations

MatrixEquations.gsylv — Function
X = gsylv(A,B,C,D,E)

Solve the generalized Sylvester matrix equation

          AXB + CXD = E

using a generalized Schur form based approach. A, B, C and D are square matrices. The pencils A-λC and D+λB must be regular and must not have common eigenvalues.

The following particular cases are also adressed:

X = gsylv(A,B,E)

Solve the generalized Sylvester matrix equation AXB = E.

X = gsylv(A,B,γ*I,E)  or  X = gsylv(A,B,γ,E)

Solve the generalized Sylvester matrix equation AXB +γX = E.

X = gsylv(A,B,γ*I,D,E)  or  X = gsylv(A,B,γ,D,E)

Solve the generalized Sylvester matrix equation AXB +γXD = E.

X = gsylv(A,B,C,δ*I,E)  or  X = gsylv(A,B,C,δ,E)

Solve the generalized Sylvester matrix equation AXB +CXδ = E.

Example

julia> A = [3. 4.; 5. 6.]
+ -3.9968e-15   -5.55112e-15
source
MatrixEquations.sylvds! — Function
sylvds!(A,B,C; adjA = false, adjB = false)

Solve the discrete Sylvester matrix equation

            op(A)Xop(B) + X =  C,

where op(A) = A or op(A) = A' if adjA = false or adjA = true, respectively, and op(B) = B or op(B) = B' if adjB = false or adjB = true, respectively. A and B are square matrices in Schur forms, and A and -B must not have common reciprocal eigenvalues. C contains on output the solution X.

source

Generalized Sylvester Matrix Equations

MatrixEquations.gsylv — Function
X = gsylv(A,B,C,D,E)

Solve the generalized Sylvester matrix equation

          AXB + CXD = E

using a generalized Schur form based approach. A, B, C and D are square matrices. The pencils A-λC and D+λB must be regular and must not have common eigenvalues.

The following particular cases are also adressed:

X = gsylv(A,B,E)

Solve the generalized Sylvester matrix equation AXB = E.

X = gsylv(A,B,γ*I,E)  or  X = gsylv(A,B,γ,E)

Solve the generalized Sylvester matrix equation AXB +γX = E.

X = gsylv(A,B,γ*I,D,E)  or  X = gsylv(A,B,γ,D,E)

Solve the generalized Sylvester matrix equation AXB +γXD = E.

X = gsylv(A,B,C,δ*I,E)  or  X = gsylv(A,B,C,δ,E)

Solve the generalized Sylvester matrix equation AXB +CXδ = E.

Example

julia> A = [3. 4.; 5. 6.]
 2×2 Array{Float64,2}:
  3.0  4.0
  5.0  6.0
@@ -78,7 +78,7 @@
 julia> A*X*B + C*X*D - E
 2×2 Array{Float64,2}:
  4.44089e-16  8.88178e-16
- 6.66134e-16  0.0
source
MatrixEquations.gsylvs! — Function
X = gsylvs!(A,B,C,D,E; adjAC=false, adjBD=false, CASchur = false, DBSchur = false)

Solve the generalized Sylvester matrix equation

            op1(A)Xop2(B) + op1(C)Xop2(D) = E,

where A, B, C and D are square matrices, and

op1(A) = A and op1(C) = C if adjAC = false;

op1(A) = A' and op1(C) = C' if adjAC = true;

op2(B) = B and op2(D) = D if adjBD = false;

op2(B) = B' and op2(D) = D' if adjBD = true.

The matrix pair (A,C) is in a generalized real or complex Schur form. The matrix pair (B,D) is in a generalized real or complex Schur form if DBSchur = false or the matrix pair (D,B) is in a generalized real or complex Schur form if DBSchur = true. The pencils A-λC and D+λB must be regular and must not have common eigenvalues.

source

Sylvester Systems of Matrix Equation

MatrixEquations.sylvsys — Function
(X,Y) = sylvsys(A,B,C,D,E,F)

Solve the Sylvester system of matrix equations

            AX + YB = C
+ 6.66134e-16  0.0
source
MatrixEquations.gsylvs! — Function
X = gsylvs!(A,B,C,D,E; adjAC=false, adjBD=false, CASchur = false, DBSchur = false)

Solve the generalized Sylvester matrix equation

            op1(A)Xop2(B) + op1(C)Xop2(D) = E,

where A, B, C and D are square matrices, and

op1(A) = A and op1(C) = C if adjAC = false;

op1(A) = A' and op1(C) = C' if adjAC = true;

op2(B) = B and op2(D) = D if adjBD = false;

op2(B) = B' and op2(D) = D' if adjBD = true.

The matrix pair (A,C) is in a generalized real or complex Schur form. The matrix pair (B,D) is in a generalized real or complex Schur form if DBSchur = false or the matrix pair (D,B) is in a generalized real or complex Schur form if DBSchur = true. The pencils A-λC and D+λB must be regular and must not have common eigenvalues.

source

Sylvester Systems of Matrix Equation

MatrixEquations.sylvsys — Function
(X,Y) = sylvsys(A,B,C,D,E,F)

Solve the Sylvester system of matrix equations

            AX + YB = C
             DX + YE = F,

where (A,D), (B,E) are pairs of square matrices of the same size. The pencils A-λD and -B+λE must be regular and must not have common eigenvalues.

Example

julia> A = [3. 4.; 5. 6.]
 2×2 Array{Float64,2}:
  3.0  4.0
@@ -129,8 +129,8 @@
 julia> D*X + Y*E - F
 2×2 Array{Float64,2}:
   1.33227e-15  -2.22045e-15
- -4.44089e-16   4.44089e-16
source
MatrixEquations.sylvsyss! — Function
(X,Y) = sylvsyss!(A,B,C,D,E,F)

Solve the Sylvester system of matrix equations

            AX + YB = C
-            DX + YE = F,

where (A,D), (B,E) are pairs of square matrices of the same size in generalized Schur forms. The pencils A-λD and -B+λE must be regular and must not have common eigenvalues. The computed solution (X,Y) is contained in (C,F).

Note: This is an enhanced interface to the LAPACK.tgsyl! function to also cover the case when A, B, D and E are real matrices and C and F are complex matrices.

source
MatrixEquations.dsylvsys — Function
(X,Y) = dsylvsys(A,B,C,D,E,F)

Solve the dual Sylvester system of matrix equations

   AX + DY = C
+ -4.44089e-16   4.44089e-16
source
MatrixEquations.sylvsyss! — Function
(X,Y) = sylvsyss!(A,B,C,D,E,F)

Solve the Sylvester system of matrix equations

            AX + YB = C
+            DX + YE = F,

where (A,D), (B,E) are pairs of square matrices of the same size in generalized Schur forms. The pencils A-λD and -B+λE must be regular and must not have common eigenvalues. The computed solution (X,Y) is contained in (C,F).

Note: This is an enhanced interface to the LAPACK.tgsyl! function to also cover the case when A, B, D and E are real matrices and C and F are complex matrices.

source
MatrixEquations.dsylvsys — Function
(X,Y) = dsylvsys(A,B,C,D,E,F)

Solve the dual Sylvester system of matrix equations

   AX + DY = C
    XB + YE = F ,

where (A,D), (B,E) are pairs of square matrices of the same size. The pencils A-λD and -B+λE must be regular and must not have common eigenvalues.

Example

julia> A = [3. 4.; 5. 6.]
 2×2 Array{Float64,2}:
  3.0  4.0
@@ -181,5 +181,5 @@
 julia> X*B + Y*E - F
 2×2 Array{Float64,2}:
  -8.88178e-16   0.0
-  8.88178e-16  -4.44089e-16
source
MatrixEquations.dsylvsyss! — Function
(X,Y) = dsylvsyss!(A,B,C,D,E,F)

Solve the dual Sylvester system of matrix equations

A'X + D'Y = C
-XB' + YE' = F,

where (A,D), (B,E) are pairs of square matrices of the same size in generalized Schur forms. The pencils A-λD and -B+λE must be regular and must not have common eigenvalues. The computed solution (X,Y) is contained in (C,F).

source
+ 8.88178e-16 -4.44089e-16source
MatrixEquations.dsylvsyss! — Function
(X,Y) = dsylvsyss!(A,B,C,D,E,F)

Solve the dual Sylvester system of matrix equations

A'X + D'Y = C
+XB' + YE' = F,

where (A,D), (B,E) are pairs of square matrices of the same size in generalized Schur forms. The pencils A-λD and -B+λE must be regular and must not have common eigenvalues. The computed solution (X,Y) is contained in (C,F).

source
diff --git a/dev/sylvkr.html b/dev/sylvkr.html index 6044d38..3da980e 100644 --- a/dev/sylvkr.html +++ b/dev/sylvkr.html @@ -1,4 +1,4 @@ -Matrix Equation Solvers using Kronecker-product Expansions · MatrixEquations.jl
GitHub

Matrix Equation Solvers using Kronecker-product Expansions

MatrixEquations.sylvckr — Function
X = sylvckr(A,B,C)

Solve the continuous Sylvester matrix equation

            AX + XB = C

using the Kronecker product expansion of equations. A and B are square matrices, and A and -B must not have common eigenvalues. This function is not recommended for large order matrices.

source
MatrixEquations.sylvdkr — Function
X = sylvdkr(A,B,C)

Solve the discrete Sylvester matrix equation

            AXB + X = C

using the Kronecker product expansion of equations. A and B are square matrices, and A and -B must not have common reciprocal eigenvalues. This function is not recommended for large order matrices.

source
MatrixEquations.gsylvkr — Function
X = gsylvkr(A,B,C,D,E)

Solve the generalized Sylvester matrix equation

            AXB + CXD = E

using the Kronecker product expansion of equations. A, B, C and D are square matrices. The pencils A-λC and D+λB must be regular and must not have common eigenvalues. This function is not recommended for large order matrices.

source
MatrixEquations.sylvsyskr — Function
sylvsyskr(A,B,C,D,E,F) -> (X,Y)

Solve the Sylvester system of matrix equations

            AX + YB = C
-            DX + YE = F

using the Kronecker product expansion of equations. (A,D), (B,E) are pairs of square matrices of the same size. The pencils A-λD and -B+λE must be regular and must not have common eigenvalues. This function is not recommended for large order matrices.

source
MatrixEquations.dsylvsyskr — Function
dsylvsyskr(A,B,C,D,E,F) -> (X,Y)

Solve the dual Sylvester system of matrix equations

   AX + DY = C
-   XB + YE = F

using the Kronecker product expansion of equations. (A,D), (B,E) are pairs of square matrices of the same size. The pencils A-λD and -B+λE must be regular and must not have common eigenvalues. This function is not recommended for large order matrices.

source
MatrixEquations.tlyapckr — Function
X = tlyapckr(A,C, isig = 1; atol::Real=0, rtol::Real=atol>0 ? 0 : N*ϵ)

Compute for isig = ±1 a solution of the the continuous T-Lyapunov matrix equation

            A*X + isig*transpose(X)*transpose(A) + C = 0

using the Kronecker product expansion of equations. A and C are m×n and m×m matrices, respectively, and X is an n×m matrix. The matrix C must be symmetric if isig = 1 and skew-symmetric if isig = -1. atol and rtol are the absolute and relative tolerances, respectively, used for rank computation. The default relative tolerance is N*ϵ, where N = 4*min(m,n)^2 and ϵ is the machine precision of the element type of A. This function is not recommended for large order matrices.

source
MatrixEquations.hlyapckr — Function
X = hlyapckr(A,C, isig = 1; atol::Real=0, rtol::Real=atol>0 ? 0 : N*ϵ)

Compute for isig = ±1 a solution of the continuous H-Lyapunov matrix equation

            A*X + isig*adjoint(X)*adjoint(A) + C = 0

using the Kronecker product expansion of equations. A and C are m×n and m×m matrices, respectively, and X is an n×m matrix. The matrix C must be hermitian if isig = 1 and skew-hermitian if isig = -1. atol and rtol are the absolute and relative tolerances, respectively, used for rank computation. The default relative tolerance is N*ϵ, where N = 4*min(m,n)^2 and ϵ is the machine precision of the element type of A. This function is not recommended for large order matrices.

source
MatrixEquations.tsylvckr — Function
X = tsylvckr(A,B,C; atol::Real=0, rtol::Real=atol>0 ? 0 : N*ϵ)

Compute a solution of the continuous T-Sylvester matrix equation

            A*X + transpose(X)*B = C

using the Kronecker product expansion of equations. A, B and C are m×n, n×m and m×m matrices, respectively, and X is an n×m matrix. atol and rtol are the absolute and relative tolerances, respectively, used for rank computation. The default relative tolerance is N*ϵ, where N = 4*min(m,n)^2 and ϵ is the machine precision of the element type of A. This function is not recommended for large order matrices.

source
MatrixEquations.hsylvckr — Function
X = hsylvckr(A,B,C; atol::Real=0, rtol::Real=atol>0 ? 0 : N*ϵ)

Compute a solution of the continuous H-Sylvester matrix equation

            A*X + adjoint(X)*B = C

using the Kronecker product expansion of equations. A, B and C are m×n, n×m and m×m matrices, respectively, and X is an n×m matrix. atol and rtol are the absolute and relative tolerances, respectively, used for rank computation. The default relative tolerance is N*ϵ, where N = 4*min(m,n)^2 and ϵ is the machine precision of the element type of A. This function is not recommended for large order matrices.

source
MatrixEquations.csylvckr — Function
X = csylvckr(A,B,C)

Solve the continuous C-Sylvester matrix equation

            A*X + conj(X)*B = C

using the Kronecker product expansion of equations. A, B and C are m×m, n×n and m×n matrices, respectively, and X is an m×n matrix. This function is not recommended for large order matrices.

source
MatrixEquations.tsylvdkr — Function
X = tsylvdkr(A,B,C)

Solve the discrete T-Sylvester matrix equation

            A*transpose(X)*B + X = C

using the Kronecker product expansion of equations. A, B and C are m×n matrices and X is an m×n matrix. This function is not recommended for large order matrices.

source
MatrixEquations.hsylvdkr — Function
X = hsylvdkr(A,B,C)

Solve the discrete H-Sylvester matrix equation

            A*adjoint(X)*B + X = C

using the Kronecker product expansion of equations. A, B and C are m×n matrices and X is an m×n matrix. This function is not recommended for large order matrices.

source
MatrixEquations.csylvdkr — Function
X = csylvdkr(A,B,C)

Solve the discrete C-Sylvester matrix equation

            A*conj(X)*B + X = C

using the Kronecker product expansion of equations. A, B and C are m×m, n×n and m×n matrices, respectively, and X is an m×n matrix. This function is not recommended for large order matrices.

source
+Matrix Equation Solvers using Kronecker-product Expansions · MatrixEquations.jl
GitHub

Matrix Equation Solvers using Kronecker-product Expansions

MatrixEquations.sylvckr — Function
X = sylvckr(A,B,C)

Solve the continuous Sylvester matrix equation

            AX + XB = C

using the Kronecker product expansion of equations. A and B are square matrices, and A and -B must not have common eigenvalues. This function is not recommended for large order matrices.

source
MatrixEquations.sylvdkr — Function
X = sylvdkr(A,B,C)

Solve the discrete Sylvester matrix equation

            AXB + X = C

using the Kronecker product expansion of equations. A and B are square matrices, and A and -B must not have common reciprocal eigenvalues. This function is not recommended for large order matrices.

source
MatrixEquations.gsylvkr — Function
X = gsylvkr(A,B,C,D,E)

Solve the generalized Sylvester matrix equation

            AXB + CXD = E

using the Kronecker product expansion of equations. A, B, C and D are square matrices. The pencils A-λC and D+λB must be regular and must not have common eigenvalues. This function is not recommended for large order matrices.

source
MatrixEquations.sylvsyskr — Function
sylvsyskr(A,B,C,D,E,F) -> (X,Y)

Solve the Sylvester system of matrix equations

            AX + YB = C
+            DX + YE = F

using the Kronecker product expansion of equations. (A,D), (B,E) are pairs of square matrices of the same size. The pencils A-λD and -B+λE must be regular and must not have common eigenvalues. This function is not recommended for large order matrices.

source
MatrixEquations.dsylvsyskr — Function
dsylvsyskr(A,B,C,D,E,F) -> (X,Y)

Solve the dual Sylvester system of matrix equations

   AX + DY = C
+   XB + YE = F

using the Kronecker product expansion of equations. (A,D), (B,E) are pairs of square matrices of the same size. The pencils A-λD and -B+λE must be regular and must not have common eigenvalues. This function is not recommended for large order matrices.

source
MatrixEquations.tlyapckr — Function
X = tlyapckr(A,C, isig = 1; atol::Real=0, rtol::Real=atol>0 ? 0 : N*ϵ)

Compute for isig = ±1 a solution of the the continuous T-Lyapunov matrix equation

            A*X + isig*transpose(X)*transpose(A) + C = 0

using the Kronecker product expansion of equations. A and C are m×n and m×m matrices, respectively, and X is an n×m matrix. The matrix C must be symmetric if isig = 1 and skew-symmetric if isig = -1. atol and rtol are the absolute and relative tolerances, respectively, used for rank computation. The default relative tolerance is N*ϵ, where N = 4*min(m,n)^2 and ϵ is the machine precision of the element type of A. This function is not recommended for large order matrices.

source
MatrixEquations.hlyapckr — Function
X = hlyapckr(A,C, isig = 1; atol::Real=0, rtol::Real=atol>0 ? 0 : N*ϵ)

Compute for isig = ±1 a solution of the continuous H-Lyapunov matrix equation

            A*X + isig*adjoint(X)*adjoint(A) + C = 0

using the Kronecker product expansion of equations. A and C are m×n and m×m matrices, respectively, and X is an n×m matrix. The matrix C must be hermitian if isig = 1 and skew-hermitian if isig = -1. atol and rtol are the absolute and relative tolerances, respectively, used for rank computation. The default relative tolerance is N*ϵ, where N = 4*min(m,n)^2 and ϵ is the machine precision of the element type of A. This function is not recommended for large order matrices.

source
MatrixEquations.tsylvckr — Function
X = tsylvckr(A,B,C; atol::Real=0, rtol::Real=atol>0 ? 0 : N*ϵ)

Compute a solution of the continuous T-Sylvester matrix equation

            A*X + transpose(X)*B = C

using the Kronecker product expansion of equations. A, B and C are m×n, n×m and m×m matrices, respectively, and X is an n×m matrix. atol and rtol are the absolute and relative tolerances, respectively, used for rank computation. The default relative tolerance is N*ϵ, where N = 4*min(m,n)^2 and ϵ is the machine precision of the element type of A. This function is not recommended for large order matrices.

source
MatrixEquations.hsylvckr — Function
X = hsylvckr(A,B,C; atol::Real=0, rtol::Real=atol>0 ? 0 : N*ϵ)

Compute a solution of the continuous H-Sylvester matrix equation

            A*X + adjoint(X)*B = C

using the Kronecker product expansion of equations. A, B and C are m×n, n×m and m×m matrices, respectively, and X is an n×m matrix. atol and rtol are the absolute and relative tolerances, respectively, used for rank computation. The default relative tolerance is N*ϵ, where N = 4*min(m,n)^2 and ϵ is the machine precision of the element type of A. This function is not recommended for large order matrices.

source
MatrixEquations.csylvckr — Function
X = csylvckr(A,B,C)

Solve the continuous C-Sylvester matrix equation

            A*X + conj(X)*B = C

using the Kronecker product expansion of equations. A, B and C are m×m, n×n and m×n matrices, respectively, and X is an m×n matrix. This function is not recommended for large order matrices.

source
MatrixEquations.tsylvdkr — Function
X = tsylvdkr(A,B,C)

Solve the discrete T-Sylvester matrix equation

            A*transpose(X)*B + X = C

using the Kronecker product expansion of equations. A, B and C are m×n matrices and X is an m×n matrix. This function is not recommended for large order matrices.

source
MatrixEquations.hsylvdkr — Function
X = hsylvdkr(A,B,C)

Solve the discrete H-Sylvester matrix equation

            A*adjoint(X)*B + X = C

using the Kronecker product expansion of equations. A, B and C are m×n matrices and X is an m×n matrix. This function is not recommended for large order matrices.

source
MatrixEquations.csylvdkr — Function
X = csylvdkr(A,B,C)

Solve the discrete C-Sylvester matrix equation

            A*conj(X)*B + X = C

using the Kronecker product expansion of equations. A, B and C are m×m, n×n and m×n matrices, respectively, and X is an m×n matrix. This function is not recommended for large order matrices.

source