Verify LDL-Decomposition of A

Pedigree

The pedigree used during the lecture corresponds to

suppressPackageStartupMessages( library(pedigreemm) )
n_nr_ani_ped <- 5
n_nr_parent <- 3
tbl_ped_simple <- dplyr::data_frame(Calf = c(1:n_nr_ani_ped),
                             Sire = c(NA, NA, NA, 1, 3),
                             Dam  = c(NA, NA, NA, 2, 2))
### # pedigreemm
(ped_simple <- pedigree(sire = tbl_ped_simple$Sire, dam = tbl_ped_simple$Dam, label = as.character(1:n_nr_ani_ped)))

Numerator Relationship Matrix

(matA_simple <- as.matrix(getA(ped = ped_simple)))

    1   2   3    4    5
1 1.0 0.0 0.0 0.50 0.00
2 0.0 1.0 0.0 0.50 0.50
3 0.0 0.0 1.0 0.00 0.50
4 0.5 0.5 0.0 1.00 0.25
5 0.0 0.5 0.5 0.25 1.00

Matrices L and D

(matD <- diag(Dmat(ped = ped_simple), n_nr_ani_ped))

     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0  0.0  0.0
[2,]    0    1    0  0.0  0.0
[3,]    0    0    1  0.0  0.0
[4,]    0    0    0  0.5  0.0
[5,]    0    0    0  0.0  0.5

### # get diagonal matrix S with diagonal elements being sqrt of elements of matD
matS <- sqrt(matD)
### # LDL decomposition based on cholesky
matR <- t(chol(matA_simple));matR

    1   2   3         4         5
1 1.0 0.0 0.0 0.0000000 0.0000000
2 0.0 1.0 0.0 0.0000000 0.0000000
3 0.0 0.0 1.0 0.0000000 0.0000000
4 0.5 0.5 0.0 0.7071068 0.0000000
5 0.0 0.5 0.5 0.0000000 0.7071068

(matL <- matR %*% solve(matS))

  [,1] [,2] [,3] [,4] [,5]
1  1.0  0.0  0.0    0    0
2  0.0  1.0  0.0    0    0
3  0.0  0.0  1.0    0    0
4  0.5  0.5  0.0    1    0
5  0.0  0.5  0.5    0    1

Verify

(matA_verify <- matL %*% matD %*% t(matL))

    1   2   3    4    5
1 1.0 0.0 0.0 0.50 0.00
2 0.0 1.0 0.0 0.50 0.50
3 0.0 0.0 1.0 0.00 0.50
4 0.5 0.5 0.0 1.00 0.25
5 0.0 0.5 0.5 0.25 1.00

Check

matA_simple - matA_verify

  1 2 3 4 5
1 0 0 0 0 0
2 0 0 0 0 0
3 0 0 0 0 0
4 0 0 0 0 0
5 0 0 0 0 0

LS0tCnRpdGxlOiAiVmVyaWZ5IExETC1EZWNvbXBvc2l0aW9uIG9mIEEiCmF1dGhvcjogIlBldGVyIHZvbiBSb2hyIgpkYXRlOiAgICIyMDE4LTExLTE2IgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgojIyBQZWRpZ3JlZQoKVGhlIHBlZGlncmVlIHVzZWQgZHVyaW5nIHRoZSBsZWN0dXJlIGNvcnJlc3BvbmRzIHRvCgpgYGB7cn0Kc3VwcHJlc3NQYWNrYWdlU3RhcnR1cE1lc3NhZ2VzKCBsaWJyYXJ5KHBlZGlncmVlbW0pICkKbl9ucl9hbmlfcGVkIDwtIDUKbl9ucl9wYXJlbnQgPC0gMwp0YmxfcGVkX3NpbXBsZSA8LSBkcGx5cjo6ZGF0YV9mcmFtZShDYWxmID0gYygxOm5fbnJfYW5pX3BlZCksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgU2lyZSA9IGMoTkEsIE5BLCBOQSwgMSwgMyksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgRGFtICA9IGMoTkEsIE5BLCBOQSwgMiwgMikpCiMjIyAjIHBlZGlncmVlbW0KKHBlZF9zaW1wbGUgPC0gcGVkaWdyZWUoc2lyZSA9IHRibF9wZWRfc2ltcGxlJFNpcmUsIGRhbSA9IHRibF9wZWRfc2ltcGxlJERhbSwgbGFiZWwgPSBhcy5jaGFyYWN0ZXIoMTpuX25yX2FuaV9wZWQpKSkKYGBgCgoKIyMgTnVtZXJhdG9yIFJlbGF0aW9uc2hpcCBNYXRyaXgKCmBgYHtyfQoobWF0QV9zaW1wbGUgPC0gYXMubWF0cml4KGdldEEocGVkID0gcGVkX3NpbXBsZSkpKQpgYGAKCgojIyBNYXRyaWNlcyBMIGFuZCBECgpgYGB7cn0KKG1hdEQgPC0gZGlhZyhEbWF0KHBlZCA9IHBlZF9zaW1wbGUpLCBuX25yX2FuaV9wZWQpKQpgYGAKCmBgYHtyfQojIyMgIyBnZXQgZGlhZ29uYWwgbWF0cml4IFMgd2l0aCBkaWFnb25hbCBlbGVtZW50cyBiZWluZyBzcXJ0IG9mIGVsZW1lbnRzIG9mIG1hdEQKbWF0UyA8LSBzcXJ0KG1hdEQpCiMjIyAjIExETCBkZWNvbXBvc2l0aW9uIGJhc2VkIG9uIGNob2xlc2t5Cm1hdFIgPC0gdChjaG9sKG1hdEFfc2ltcGxlKSkKKG1hdEwgPC0gbWF0UiAlKiUgc29sdmUobWF0UykpCmBgYAoKCiMjIFZlcmlmeQpgYGB7cn0KKG1hdEFfdmVyaWZ5IDwtIG1hdEwgJSolIG1hdEQgJSolIHQobWF0TCkpCmBgYAoKIyMgQ2hlY2sKYGBge3J9Cm1hdEFfc2ltcGxlIC0gbWF0QV92ZXJpZnkKYGBgCgoKCgoK