Ridge estimation for high-dimensional precision matrices with known chordal support
Calculates various ridge estimators for high-dimensional precision matrices with known support. This support should form a chordal graph. If the provided support is not chordal, the function makes it so.
ridgePchordal( S, lambda, zeros, cliques = list(), separators = list(), target = default.target(S), type = "Alt", optimizer = "nlm", grad = FALSE, verbose = TRUE, ... )
S |
Sample covariance |
lambda |
A |
zeros |
|
cliques |
A |
separators |
A |
target |
A target |
type |
A |
optimizer |
A |
grad |
A |
verbose |
A |
... |
Sister function to the ridgeP-function, incorporating a chordal
zero structure of the precision matrix.
The loss function for type="ArchII" is:
\log(| \mathbf{Ω} |) - \mbox{tr}( \mathbf{S} \mathbf{Ω} ) + λ \big\{ \log(| \mathbf{Ω} |) - \mbox{tr}[ (\mathbf{S} + (1+λ) \mathbf{I}_{p \times p}) \mathbf{Ω} ] \big\}.
For type="ArchI" it is:
(1-λ) \big[ \log(| \mathbf{Ω} |) - \mbox{tr}( \mathbf{S} \mathbf{Ω} ) \big] + λ \big[ \log(| \mathbf{Ω} |) - \mbox{tr}( \mathbf{Ω} ) \big],
which is obtained from:
\log(| \mathbf{Ω} |) - \mbox{tr}( \mathbf{S} \mathbf{Ω} ) + ν \big[ \log(| \mathbf{Ω} |) - \mbox{tr}( \mathbf{Ω} ) \big]
by division of (1+ν) and writing λ = ν / (1 + ν).
An explicit expression for the minimizer of the loss functions implied by the
archetypal ridge estimators (type="ArchI" and type="ArchII")
exists. For the simple case in which the graph decomposes into cliques
\mathcal{C}_1, \mathcal{C}_2 and separator \mathcal{S} the
estimator is:
\widehat{\mathbf{Ω}} = ≤ft( \begin{array}{lll} \, [\widehat{\mathbf{Ω}}^{({\mathcal{C}_1})}]_{\mathcal{C}_1 \setminus \mathcal{S}, \mathcal{C}_1 \setminus \mathcal{S}} & [\widehat{\mathbf{Ω}}^{({\mathcal{C}_1})}]_{\mathcal{C}_1 \setminus \mathcal{S}, \mathcal{S}} & \mathbf{0}_{|\mathcal{C}_1 \setminus \mathcal{S}| \times |\mathcal{C}_2 \setminus \mathcal{S}|} \\ \, [\widehat{\mathbf{Ω}}^{(\mathcal{C}_1)}]_{\mathcal{S}, \mathcal{C}_1 \setminus \mathcal{S}} & [\widehat{\mathbf{Ω}}^{({\mathcal{C}_1})}]_{\mathcal{S}, \mathcal{S}} + [\widehat{\mathbf{Ω}}^{({\mathcal{C}_2})}]_{\mathcal{S}, \mathcal{S}} - \widehat{\mathbf{Ω}}^{(\mathcal{S})} & [\widehat{\mathbf{Ω}}^{({\mathcal{C}_2})}]_{\mathcal{S}, \mathcal{C}_2 \setminus \mathcal{S}} \\ \, \mathbf{0}_{|\mathcal{C}_2 \setminus \mathcal{S}| \times |\mathcal{C}_1 \setminus \mathcal{S}|} & [\widehat{\mathbf{Ω}}^{({\mathcal{C}_2})}]_{\mathcal{C}_2 \setminus \mathcal{S}, \mathcal{S}} & [\widehat{\mathbf{Ω}}^{({\mathcal{C}_2})}]_{\mathcal{C}_2 \setminus \mathcal{S}, \mathcal{C}_2 \setminus \mathcal{S}} \end{array} \right),
where \widehat{\mathbf{Ω}}^{({\mathcal{C}_1})}, \widehat{\mathbf{Ω}}^{({\mathcal{C}_1})} and \widehat{\mathbf{Ω}}^{({\mathcal{S}})} are the marginal ridge ML covariance estimators for cliques \mathcal{C}_1, \mathcal{C}_2 and separator \mathcal{S}. The general form of the estimator, implemented here, is analogous to that provided in Proposition 5.9 of Lauritzen (2004). The proof that this estimator indeed optimizes the corresponding loss function is fully analogous to that of Proposition 5.6 of Lauritzen (2004).
In case, type="Alt" no explicit expression of the maximizer of the
ridge penalized log-likelihood exists. However, an initial estimator
analogous to that for type="ArchI" and type="ArchII" can be
defined. In various boundary cases (λ=0, λ=∞,
and \mathcal{S} = \emptyset) this initial estimator actually optimizes
the loss function. In general, however, it does not. Nevertheless, it
functions as well-educated guess for any Newton-like optimization method:
convergence is usually achieved quickly. The Newton-like procedure optimizes
an unconstrained problem equivalent to that of the penalized log-likelihood
with known zeros for the precision matrix (see Dahl et al., 2005 for
details).
If grad=FALSE, the function returns a regularized precision
matrix with specified chordal sparsity structure.
If grad=TRUE, a list is returned comprising of i) the
estimated precision matrix, and ii) the gradients at the initial and
at the optimal (if reached) value. The gradient is returned and it can be
checked whether it is indeed (close to) zero at the optimum.
Wessel N. van Wieringen.
Dahl, J., Roychowdhury, V., Vandenberghe, L. (2005), "Maximum likelihood estimation of Gaussian graphical models: numerical implementation and topology selection", Technical report, UCLA, 2005.
Lauritzen, S.L. (2004). Graphical Models. Oxford University Press.
Miok, V., Wilting, S.M., Van Wieringen, W.N. (2016), "Ridge estimation of the VAR(1) model and its time series chain graph from multivariate time-course omics data", Biometrical Journal, 59(1), 172-191.
# obtain some (high-dimensional) data p <- 8 n <- 100 set.seed(333) Y <- matrix(rnorm(n*p), nrow = n, ncol = p) # define zero structure S <- covML(Y) S[1:3, 6:8] <- 0 S[6:8, 1:3] <- 0 zeros <- which(S==0, arr.ind=TRUE) # obtain (triangulated) support info supportP <- support4ridgeP(nNodes=p, zeros=zeros) # estimate precision matrix with known (triangulated) support Phat <- ridgePchordal(S, 0.1, zeros=supportP$zeros, cliques=supportP$cliques, separators=supportP$separators)
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.