Introduction to mbend
In mbend: Matrix Bending

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Introduction to mbend

Description

Bending non-positive-definite (symmetric) matrices to positive-definite, using weighted and unweighted methods

The mbend package is used for bending symmetric non-positive-definite matrices to positive-definite (PD). For a matrix to be invertible, it has to be PD. The methods of Jorjani et al. (2003) and Schaeffer (2014) are used in this package. The unweighted method of Schaeffer (2014) is also extended to weighted bending, with the possibility of choosing between unweighted and weighted bending.

Application

Start with loading the package library:

library(mbend)

Consider the following non-PD covariance matrix (Jorjani et al., 2003):

V = matrix(nrow=5, ncol=5, c(
  100,  95,  80,  40,  40,
   95, 100,  95,  80,  40,
   80,  95, 100,  95,  80,
   40,  80,  95, 100,  95,
   40,  40,  80,  95, 100))

The following command can be used to bend matrix V to a PD matrix:

bend(V)

The above command is equivalent to bend(inmat=V), where inmat is the argument that takes the matrix to be bent, or the following command:

bend(V, max.iter=10000, small.positive=0.0001, method="hj")

This runs the unweighted bending method of Jorjani et al. (2003) (method="hj"), with maximum 10000 number of iterations (max.iter=10000), and eigenvalues smaller than 0.0001 are replaced with this small positive value (small.positive=0.0001). Providing the default parameters, the corresponding arguments can be omitted (e.g., bend(V)).

The output object is a list of several items, listed below:

bent : The bent matrix.
init.ev : Eigenvalues of the initial (inmat) matrix.
final.ev : Eigenvalues of the bent matrix.
min.dev : min(bent - inmat).
max.dev : max(bent - inmat).
loc.min.dev : Location (indices) of min.dev element.
loc.max.dev : Location (indices) of max.dev element.
ave.dev : Average deviation (bent - inmat) of the upper triangle elements (excluding diagonal elements for correlation matrices).
AAD : Average absolute deviation of the upper triangle elements (excluding diagonal elements for correlation matrices) of bent and inmat.
Cor : Correlation between the upper triangle elements (excluding diagonal elements for correlation matrices) of bent and inmat.
RMSD : Root of mean squared deviation of the upper triangle elements (excluding diagonal elements for correlation matrices) of bent and inmat.

There might be different precision involved with different elements of a non-PD matrix. In this case, a weighted bending is recommended. Jorjani et al. (2003) used the reciprocal of the number data points in common between pairs of variables, as weights. Considering the following matrix for the number of data points in common between variables (Jorjani et al., 2003):

W = matrix(nrow=5, ncol=5, c(
  1000,  500,   20,   50,  200,
   500, 1000,  500,    5,   50,
    20,  500, 1000,   20,   20,
    50,    5,   20, 1000,  200,
   200,   50,   20,  200, 1000))

Matrix V is bent using the following command:

bend(inmat=V, wtmat=W, reciprocal=TRUE)

Using wtmat=1/W, reciprocal=FALSE, the argument reciprocal could be omitted, because FALSE is the default parameter for the argument reciprocal. For the same reason, max.iter=10000, small.positive=0.0001, method="hj" are omitted, unless different parameters are provided to these arguments.

If there is high confidence about some elements of the non-PD matrix to remain unchanged after bending, the corresponding weights are set to zero. For example, to keep the first 2 × 2 block of V unchanged during the bending procedure:

W2 = W; W2[1:2, 1:2] = 0
bend(V, W2, reciprocal=TRUE)

For weighted bending, we get extra statistics in the output:

w_gt_0 : Number of weight elements greater than 0, in the upper triangle of wtmat (for weighted bending).
wAAD : Weighted AAD (for weighted bending).
wCor : Weighted Cor (for weighted bending).
wRMSD : Weighted RMSD (for weighted bending).

To bend V using the method of Schaeffer (2014):

bend(inmat=V, method="lrs")

The method of Schaeffer (2014) does not require the argument small.positive, and this argument is ignored. This method is originally an unweighted bending method. However, in this package, it is extended to accommodate weighted bending (i.e., a combination of Schaeffer (2014) and Jorjani et al. (2003) methods). Weighted bending of V with reciprocals of W using the method of Schaeffer (2014):

bend(V, W, reciprocal=TRUE, method="lrs")

Function bend automatically considers any matrix with all diagonal elements equal to one, as a correlation matrix. Considering the correlation matrix V2 (V converted to a correlation matrix):

V2 = cov2cor(V)
bend(V2, W, reciprocal=TRUE)

For correlation matrices, diagonal elements are not used in obtaining the statistics. Because the argument method is not provided, the default parameter "hj" is used. To do the same using the method of Schaeffer (2014):

bend(V2, W, reciprocal=TRUE, method="lrs")

Bending the same correlation matrix using the unweighted Schaeffer (2014):

bend(V2, method="lrs")

References

Jorjani, H., Klie. L., & Emanuelson, U. (2000). A simple method for weighted bending of genetic (co)variance matrices. J. Dairy Sci. 86(2): 677--679. https://doi.org/10.3168/jds.S0022-0302(03)73646-7

Schaeffer, L. R. (2014). Making covariance matrices positive definite. Available at: http://animalbiosciences.uoguelph.ca/~lrs/ELARES/PDforce.pdf