Contents Previous Next Subchapters

Solves a Symmetric Block Tri-Diagonal System of Equations
Syntax [Elam] = BlockTriDiag(BCR)
[Elam] = BlockTriDiag(BCRlevel)
See Also tridiag

Description
Solves for the sequence of matrices {e(k)} in the system of equations

     / b(1)  c(2)' 0       ... 0 \ / e(1) \   / r(1) \
     | c(2)  b(2)  c(3)'       . | | e(2) |   | r(2) |
     |             .           . | |  .   | = |  .   |
     |               .           | |  .   |   |  .   |
     \ 0   ...   0    c(N)  b(N) / \ e(N) /   \ r(N) /

where each b(k) is a symmetric positive definite n by n matrix, each c(k) is an n by n matrix, each e(k) is an n by m matrix, and each r(k) is an n by m matrix.

Assumptions
The sequence of matrices {d(k)}, generated by the algorithm, are assumed to be invertible. This is guarantee if the block tri-diagonal matrix above is positive definite (see Vassilevski's article).

B
The real, or double-precision matrix B has n * N rows and n columns. The value of b(k) is equal to the following submatrix of B:
     B.row((k-1) * n + 1, n)

C
The real, or double-precision matrix C has n * N rows and n columns. The value of c(k) is equal to the following submatrix of C:
     C.row((k-1) * n + 1, n)
Note that the first n rows of C are not used.

R
The real, or double-precision matrix R has n * N rows and m columns. The value of r(k) is equal to the following submatrix of R:
     R.row((k-1) * n + 1, m)

E
The real, or double-precision matrix E has n * N rows and m columns. The value of e(k) is equal to the following submatrix of E:
     E.row((k-1) * n + 1, m)

lam
The real or double-precision scalar lam is set to the log of the determinant of the tri-diagonal matrix in the system of equations above.

level
Is an integer scalar specifying the level of tracing to do. If level < 0, no tracing is done. This is the default value for level if it is not present. The block tri-diagonal matrix above is formed in memory and referred to a T below. This can be a very large matrix and is not formed unless level > 0.
Level Label Description
level > 1 type(T) type of the matrix T
level > 1 cond(T) condition number of the matrix T
level > 1 lam block tri-diagonal calculation of logdet(T)
level > 1 logdet(T) direct calculation of logdet(T)
level > 2 E block tri-diagonal calculation of solution of equation
level > 2 T \ R direct calculation of solution of equation
level > 3 T the large block tri-diagonal matrix
Example

clear
include BlockTriDiag.oms
#
function RandomPositive(n) begin
     F = rand(n, n)
     return F' * F
end
#
# dimensions of the problem
N  = 4
n  = 2
m  = 1
#
# dimension values
T = fill(0d0, n * N, n * N)
Y = fill(0d0, n * N, m)
B = fill(0d0, n * N, n)
C = fill(0d0, n * N, n)
R = fill(0d0, n * N, m)
#
# Generate the sequences q(k), a(k), u(k), b(k), and c(k)
# as per the conditions in Lemma 7
#
# q(0) and a(0)
qkm = RandomPositive(n)
akm = rand(n, n)
for k = 1 to N begin
     # generate random u(k), q(k), a(k), r(k)
     uk = RandomPositive(n)
     qk = RandomPositive(n)
     ak = rand(n, n)
     rk = rand(n, m)
     #
     bk = uk + inv(qkm) + ak * ( qk \ ak' )
     ck = qkm \ akm'
     #
     # Fill in block tri-diagonal system of equations 
     kn = (k - 1) * n + 1
     T.blk(kn, kn, n, n) = bk
     if k > 1 then begin
          T.blk(kn, kn - n, n, n) = ck
          T.blk(kn - n, kn, n, n) = ck'
     end
     Y.blk(kn, 1, n, m) = rk
     #
     # store b(k), c(k), and r(k) in form for BlockTriDiag
     kn = (k - 1) * n + 1
     B.row(kn, n) = bk
     C.row(kn, n) = ck
     R.row(kn, n) = rk
end
#
# solve the block Tri-Diagonal System
level    = 3
[E, lam] = BlockTriDiag(B, C, R, level)

Reference
Algorithm 5 in Bell, B.M., The marginal likelihood corresponding to a discrete Gauss-Markov process IEEE Transactions on Signal Processing, March 2000.