Contents Previous Next Subchapters Current Chapters-> kron inv trace det logdet rank cond null orth pinv cholesky lu svd qr qred eigen geneig eigsym symeig schur levinson tridiag BlockTriDiag Parent Chapters-> Omatrix6 linear BlockTriDiag Search Tools-> contents reference index search

Solves a Symmetric Block Tri-Diagonal System of Equations
 Syntax `[`E`, `lam`] = BlockTriDiag(`B`, `C`, `R`)` `[`E`, `lam`] = BlockTriDiag(`B`, `C`, `R`, `level`)` 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.