Utilities for Jordan matrices — jordan • mcompanion

Utility functions for work with the Jordan decompositions of matrices: create a block diagonal matrix of Jordan blocks, restore a matrix from its Jordan decomposition, locate specific chains.

Usage

Jordan_matrix(eigval, len.block)
from_Jordan(x, jmat, ...)
chain_ind(chainno, len.block)
chains_to_list(vectors, heights)

Arguments

eigval: eigenvalues, a numeric or complex vector.
len.block: lengths of Jordan chains, a vector of positive integers.
x: generalised eigenvectors, a matrix with one column for each (generalised) eigenvector.
jmat: a Jordan matrix.
chainno: a vector of positive integers between 1 and length(eigval) specifying which Jordan chains to locate, see Details.
...: further arguments to pass on to solve.
vectors: a matrix of generalised eigenvectors of a matrix.
heights: a vector of chain lengths, heights[i] is the length of the i-th chain.

Details

Jordan_matrix creates a Jordan matrix (block-diagonal matrix with Jordan blocks on the diagonal) whose i-th diagonal block corresponds to eigval[i] and is of size len.block[i]. If len.block is missing, Jordan_matrix returns diag(eigenvalues).

from_Jordan computes the matrix whose Jordan decomposition is represented by arguments X (chains) and J (Jordan matrix). Conceptually, the result is equivalent to \(XJX^{-1}\) but without explicitly inverting matrices (currently the result is the transpose of solve(t(x), t(x %*% jmat), ...)).

chain_ind computes the columns of specified Jordan chains in a matrix of generalised eigenvectors. It is mostly internal function. If x is a matrix whose columns are generalised eigenvectors and the i-th Jordan chain is of length len.block[i], then this function gives the column numbers of x containing the specified chains. Note that chain_ind is not able to deduce the total number of eigenvalues. It is therefore an error to omit argument len.block when calling it.

chains_to_list converts the matrix vectors into a list of matrices. The i-th element of this list is a matrix whose columns are the vectors in the i-th chain.

Value

for Jordan_matrix, a matrix with the specified Jordan blocks on its diagonal.

for from_Jordan, the matrix with the specified Jordan decomposition.

for chain_ind, a vector of positive integers giving the columns of the requested chains.

for chains_to_list, a list of matrices.

Author

Georgi N. Boshnakov

Level

Examples

## single Jordan blocks
Jordan_matrix(4, 2) 
#>      [,1] [,2]
#> [1,]    4    1
#> [2,]    0    4
Jordan_matrix(5, 3)
#>      [,1] [,2] [,3]
#> [1,]    5    1    0
#> [2,]    0    5    1
#> [3,]    0    0    5
Jordan_matrix(6, 1)
#>      [,1]
#> [1,]    6
## a matrix with the above 3 blocks
Jordan_matrix(c(4, 5, 6), c(2, 3, 1))
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,]    4    1    0    0    0    0
#> [2,]    0    4    0    0    0    0
#> [3,]    0    0    5    1    0    0
#> [4,]    0    0    0    5    1    0
#> [5,]    0    0    0    0    5    0
#> [6,]    0    0    0    0    0    6

## a matrix with a 2x2 Jordan block for eval 1 and two simple 0 eval's
m <- make_mcmatrix(eigval = c(1), co = cbind(c(1,1,1,1), c(0,1,0,0)),
                     dim = 4, len.block = c(2))
m
#>      [,1] [,2] [,3] [,4]
#> [1,]    0    1    0    0
#> [2,]   -1    2    0    0
#> [3,]    0    1    0    0
#> [4,]    0    1    0    0
m.X <- cbind(c(1,1,1,1), c(0,1,0,0), c(0,0,1,0), c(0,0,0,1))
m.X
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    0    0    0
#> [2,]    1    1    0    0
#> [3,]    1    0    1    0
#> [4,]    1    0    0    1
m.J <- cbind(c(1,0,0,0), c(1,1,0,0), rep(0,4), rep(0,4))
m.J
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    1    0    0
#> [2,]    0    1    0    0
#> [3,]    0    0    0    0
#> [4,]    0    0    0    0

from_Jordan(m.X, m.J)          # == m
#>      [,1] [,2] [,3] [,4]
#> [1,]    0    1    0    0
#> [2,]   -1    2    0    0
#> [3,]    0    1    0    0
#> [4,]    0    1    0    0
m.X %*% m.J %*% solve(m.X) # == m
#>      [,1] [,2] [,3] [,4]
#> [1,]    0    1    0    0
#> [2,]   -1    2    0    0
#> [3,]    0    1    0    0
#> [4,]    0    1    0    0
all(m == from_Jordan(m.X, m.J)) && all(m == m.X %*% m.J %*% solve(m.X))
#> [1] TRUE
## TRUE

## which column(s) in m.X correspond to 1st Jordan block?
chain_ind(1, c(2,1,1)) # c(1, 2) since 2x2 Jordan block
#> [1] 1 2
    
## which column(s) in m.X correspond to 2nd Jordan block?
chain_ind(2, c(2,1,1)) # 3, simple eval
#> [1] 3

## which column(s) in m.X correspond to 1st and 2nd Jordan blocks?
chain_ind(c(1, 2), c(2,1,1)) # c(1,2,3)
#> [1] 1 2 3
## non-contiguous subset are ok:
chain_ind(c(1, 3), c(2,1,1)) # c(1,2,4)
#> [1] 1 2 4

## split the chains into a list of matrices
chains_to_list(m.X, c(2,1,1))
#> [[1]]
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    1    1
#> [3,]    1    0
#> [4,]    1    0
#> 
#> [[2]]
#>      [,1]
#> [1,]    0
#> [2,]    0
#> [3,]    1
#> [4,]    0
#> 
#> [[3]]
#>      [,1]
#> [1,]    0
#> [2,]    0
#> [3,]    0
#> [4,]    1
#> 
    
m.X %*% m.J
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    1    0    0
#> [2,]    1    2    0    0
#> [3,]    1    1    0    0
#> [4,]    1    1    0    0
m %*% m.X    # same
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    1    0    0
#> [2,]    1    2    0    0
#> [3,]    1    1    0    0
#> [4,]    1    1    0    0
all(m.X %*% m.J == m %*% m.X)    # TRUE
#> [1] TRUE
        
m %*% c(1,1,1,1)    # = c(1,1,1,1),  evec for eigenvalue 1
#>      [,1]
#> [1,]    1
#> [2,]    1
#> [3,]    1
#> [4,]    1
m %*% c(0,1,0,0)    # gen.e.v. for eigenvalue 1
#>      [,1]
#> [1,]    1
#> [2,]    2
#> [3,]    1
#> [4,]    1
## indeed:
all( m %*% c(0,1,0,0) == c(0,1,0,0) + c(1,1,1,1) ) # TRUE    
#> [1] TRUE

## m X = X jordan.block
cbind(c(1,1,1,1), c(0,1,0,0)) %*% cbind(c(1,0), c(1,1))
#>      [,1] [,2]
#> [1,]    1    1
#> [2,]    1    2
#> [3,]    1    1
#> [4,]    1    1
m %*% cbind(c(1,1,1,1), c(0,1,0,0))
#>      [,1] [,2]
#> [1,]    1    1
#> [2,]    1    2
#> [3,]    1    1
#> [4,]    1    1