Solve Rational Matrix Equation

Aditi Tiwari

Given a symmetrix positive definite matrix Q and a non-singular matrix L, Find symmetric positive definite solution X such that X = Q + L (X inv) L^T.

Solve Rational Matrix

library(SolveRationalMatrixEquation)

Q <- rbind(c(2,-1,0),c(-1,2,1),c(0,1,2))
L <- rbind(c(2,-2,18),c(2,1,0),c(1,2,0))

##QR Decomposition
QRdecompose(L)
## $Q
##            [,1]       [,2]       [,3]
## [1,] -0.6666667  0.6666667 -0.3333333
## [2,] -0.6666667 -0.3333333  0.6666667
## [3,] -0.3333333 -0.6666667 -0.6666667
## 
## $R
##               [,1]          [,2] [,3]
## [1,] -3.000000e+00  2.220446e-16  -12
## [2,] -1.165734e-16 -3.000000e+00   12
## [3,]  1.554312e-16  0.000000e+00   -6
##LQ Decompostion
LQdecompose(L)
## $L
##             [,1]          [,2]          [,3]
## [1,] -18.2208672 -2.532680e-16  4.270531e-16
## [2,]  -0.1097643 -2.233372e+00  0.000000e+00
## [3,]   0.1097643 -1.796408e+00 -1.326978e+00
## 
## $Q
##            [,1]       [,2]        [,3]
## [1,] -0.1097643  0.1097643 -0.98787834
## [2,] -0.8901121 -0.4531480  0.04855157
## [3,]  0.4423259 -0.8846517 -0.14744196
##Solve Rational Matrix
X <- sol.rationalmatrix.euqation(Q, L)

#Checking the results
X - (Q + L %*% solve(X) %*% t(L))
##           [,1]        [,2]          [,3]
## [1,] 0.5000000  0.13397460 -1.776357e-15
## [2,] 0.1339746  0.16666667 -5.719096e-02
## [3,] 0.0000000 -0.05719096 -6.666667e-01