Cholesky Decomposition

Cholesky Decomposition for Sampling Correlated Parameters in Probabilistic Sensitivity Analysis

In a health economic model, the Cholesky decomposition is the usual way to make sampled parameters move together in the way the evidence says they should. It splits the parameters' covariance matrix into a lower-triangular factor, which turns independent random draws into draws that carry the variances and correlations of a fitted regression, a parametric survival curve or a set of related inputs, so that each iteration of a probabilistic sensitivity analysis uses a coherent set of values. The Canadian economic evaluation guidelines require it when regression coefficients feed a model, and the Dutch guideline names it for the parameters of survival curves fitted to Kaplan-Meier data.

Why correlated parameters need joint sampling

Many model inputs are estimated together and are therefore uncertain together. The ISPOR-SMDM report on parameter estimation and uncertainty (Task Force-6) notes that when regression captures the effect of patient characteristics, uncertainty in the derived parameter is defined by uncertainty and correlation in the coefficients, and that multivariate normality is appropriate for the regression's linear predictor. Its best-practice statement VI-10 says that correlation among parameters should be considered, and that jointly estimated parameters such as regression coefficients carry direct evidence on correlation which should be reflected in the analysis.

Survival extrapolation raises the same issue. NICE DSU Technical Support Document 14 states that, when patient-level data are available, parameter uncertainty can be taken into account using the variance-covariance matrices of the fitted parametric models. The parameters of a curve such as a Weibull model are estimated from the same data, so they are uncertain together. Neine and Curran give a third case. In a herpes zoster vaccination model stratified by age, sampling each age group's incidence independently can produce iterations in which an older group has lower incidence than a younger one, contrary to the epidemiology. They also report that a review of cost-effectiveness models conducted for NICE found that most ignored potential correlation between inputs.

The factorisation and when it exists

Higham defines a symmetric matrix as positive definite when the quadratic form $x^{\mathsf T}Ax$ is positive for every non-zero vector $x$, or equivalently when all its eigenvalues are positive. Such a matrix can be written as a triangular factor times its own transpose. With a positive diagonal the factor is unique, as follows.

$$A = L L^{\mathsf T}$$

where $A$ is an $n \times n$ symmetric positive-definite matrix, $L$ is lower triangular with positive diagonal entries and $L^{\mathsf T}$ is its transpose. Higham also writes it as $A = R^{\mathsf T}R$ with $R = L^{\mathsf T}$ upper triangular. A covariance matrix of model parameters is symmetric, and when no parameter is an exact linear combination of the others it is positive definite.

The entries are found one column at a time, each from entries already computed. Equating elements of $A$ and $LL^{\mathsf T}$ gives the following for column $j$ and rows $i > j$.

$$\ell_{jj} = \sqrt{ a_{jj} - \sum_{k=1}^{j-1} \ell_{jk}^{2} }, \qquad \ell_{ij} = \frac{1}{\ell_{jj}} \left( a_{ij} - \sum_{k=1}^{j-1} \ell_{ik}\ell_{jk} \right)$$

where $a_{ij}$ is the element of $A$ in row $i$ and column $j$, $\ell_{ij}$ is the corresponding element of $L$, and $k$ runs over the columns already completed. Higham notes that positive definiteness guarantees the quantity under each square root is positive, and that the algorithm needs about $n^3/3$ floating-point operations. Neine and Curran found that adding Cholesky-based correlation to their Excel model added very little computational overhead.

From independent draws to correlated parameters

The sampling step is a linear transformation that needs only the point estimates and their covariance matrix. Each iteration draws independent standard normal values and multiplies them by the Cholesky factor, as follows.

$$\theta = \mu + L z, \qquad z \sim N(0, I), \qquad \text{Var}(\theta) = L \, I \, L^{\mathsf T} = \Sigma$$

where $\theta$ is the vector of sampled parameters for one iteration, $\mu$ is the vector of point estimates, $\Sigma = LL^{\mathsf T}$ is their covariance matrix, $z$ is a vector of independent standard normal draws, $I$ is the identity matrix and $\text{Var}(\theta)$ is the covariance of the sampled vector. Because the transformation is linear, $\theta$ follows a multivariate normal distribution with mean $\mu$ and covariance $\Sigma$, the form Task Force-6 calls appropriate for a regression's linear predictor.

In practice the covariance matrix is factorised once, a new $z$ is drawn in every iteration of the Monte Carlo simulation, and every quantity the coefficients feed, such as a subgroup probability or a survival curve, is recomputed from $\theta$ before the model runs.

Worked example with two correlated logistic regression coefficients

All numbers in this example are illustrative. A logistic regression predicts the annual probability of hospital admission from the number of previous admissions. The intercept is −2.0 with standard error 0.20, the coefficient on previous admissions is 0.5 with standard error 0.08, and their estimated correlation is −0.6, so the covariance is −0.6 × 0.20 × 0.08 = −0.0096.

1. Factorise the covariance matrix. The variances are 0.20 squared = 0.04 and 0.08 squared = 0.0064. The first column gives a diagonal entry of 0.2, the square root of 0.04, and a lower entry of −0.0096 / 0.2 = −0.048. The second diagonal entry is the square root of 0.0064 − 0.002304 = 0.004096, which is 0.064, where 0.002304 is 0.048 squared.

$$L = \begin{pmatrix} 0.2 & 0 \\ -0.048 & 0.064 \end{pmatrix}, \qquad LL^{\mathsf T} = \begin{pmatrix} 0.04 & -0.0096 \\ -0.0096 & 0.0064 \end{pmatrix}$$

where $L$ is the lower-triangular Cholesky factor, and multiplying it by its transpose recovers the covariance matrix, since 0.2 × −0.048 = −0.0096 and 0.002304 + 0.004096 = 0.0064.

2. Build one draw. Suppose an iteration draws $z = (1.0, -0.5)$. The sampled intercept is −2.0 + 0.2 × 1.0 = −1.8 and the sampled coefficient is 0.5 + (−0.048 × 1.0) + (0.064 × −0.5) = 0.42. The shared first draw pushes the coefficient down when the intercept goes up, which is the negative correlation the regression estimated.

3. Compute the model input. For a patient with one previous admission the linear predictor at the point estimates is −2.0 + 0.5 = −1.5, a probability of 1 / (1 + e^1.5) = 0.182. In this iteration it is −1.8 + 0.42 = −1.38, a probability of 1 / (1 + e^1.38) = 0.201.

4. See what independence would do. The variance of this patient's linear predictor is 0.04 + 0.0064 + 2 × (−0.0096) = 0.0272, a standard error of 0.165. Sampling the coefficients independently drops the covariance term and gives 0.04 + 0.0064 = 0.0464, a standard error of 0.215, about 31% too large. The approximate 95% range for the admission probability widens from 0.139 to 0.236 with correlation to 0.128 to 0.254 without it.

5. Interpretation. Ignoring this negative correlation overstates uncertainty in the model input and spreads the cloud of results on the cost-effectiveness plane. With a positive correlation, or for a quantity that depends on a difference between parameters, independence can understate uncertainty instead. Simulating many draws and comparing the sample covariance matrix with the target confirms that the factor is being applied correctly.

Beta, gamma and other non-normal parameters

Most model inputs are not normal. Probabilities usually take a beta distribution and costs a gamma distribution. Neine and Curran describe an algorithm that uses the Cholesky factor to correlate standard normal draws, converts each correlated value to a probability with the normal cumulative distribution function, and then applies the inverse cumulative distribution function of the chosen beta, gamma or log-normal distribution. They implemented it in Excel with a Visual Basic for Applications macro.

The correlations achieved after transformation were close to, but not exactly, those specified. For targets of 0.25, 0.5, 0.75 and 0.9, the median correlations obtained with gamma marginals were 0.247, 0.496, 0.748 and 0.899. When only incidence was varied, 90% of iterations gave ICERs below USD 25,000, USD 33,000 and USD 38,000 per QALY with correlations of 0, 0.5 and 0.9, so the assumed correlation changed the decision uncertainty. For probabilities that must sum to one, such as a row of transition probabilities, the CADTH guidelines name the Dirichlet distribution, which handles the constraint without a Cholesky factor.

The orientation trap in software

Software differs in which triangle it returns. The R function chol returns the upper-triangular factor $R$ such that $R^{\mathsf T}R$ equals the input, and its documentation warns that the code does not check for symmetry. The NumPy function numpy.linalg.cholesky returns the lower-triangular factor by default and the upper one only on request.

In the worked example, using the upper factor $R = L^{\mathsf T}$ in place of $L$ would give the sampled vector covariance $RR^{\mathsf T}$, with a variance of 0.042304 for the intercept, 0.004096 for the coefficient and a covariance of −0.003072, against targets of 0.04, 0.0064 and −0.0096. The upper factor is correct when a row vector of draws multiplies it from the left, $z^{\mathsf T}R$, the form Neine and Curran use. No error appears, but a validation run comparing sample and target covariances catches it.

What a failed factorisation says about the inputs

Failure is informative. Higham describes running the algorithm as a fast test of positive definiteness, in which the matrix is declared positive definite if no zero or negative pivot appears. The LAPACK routine DPOTRF, which R calls, reports when the leading minor of a given order is not positive definite, and R's chol stops with an error for a matrix that is only positive semidefinite unless pivoting is requested.

In economic models the usual cause is a correlation matrix assembled pair by pair. Higham defines a correlation matrix as symmetric positive semidefinite with unit diagonal, and separate pairwise judgements need not meet that condition. As an illustration, three parameters with correlations of 0.9 between the first and second, 0.9 between the first and third, and 0 between the second and third have a determinant of 1 − 0.81 − 0.81 = −0.62. The factorisation gives a second pivot of 0.436, the square root of 0.19, and a lower entry of −0.81 / 0.436 = −1.858, and the third pivot would need the square root of 1 − 0.81 − 3.453 = −3.263, so it fails. Two parameters that each move closely with a third cannot be uncorrelated with each other.

Higham's nearest correlation matrix method repairs a matrix built from inconsistent data, although revisiting the elicitation or evidence behind the correlations is often the better response. Bayesian evidence synthesis meets the same constraint, and NICE DSU Technical Support Document 20 warns that independent priors on the elements of a between-studies covariance matrix can produce sampled matrices that are not positive semidefinite, and describes parameterising the correlation matrix through its Cholesky factor in multivariate meta-analysis.

Alternatives and limits of the Cholesky approach

The Cholesky approach is one of several ways to respect correlation. Task Force-6 notes that re-parameterising a model so that inputs are plausibly independent is often practical and sufficient. A baseline progression probability and a relative risk reduction can be sampled independently, with the treated probability derived as their product. When parameters come from Markov chain Monte Carlo, the posterior samples already carry their joint distribution and no factorisation is needed.

The method has limits that are easy to overlook. Neine and Curran note that the correlation is treated as fixed in every iteration, so within their method uncertainty about the correlation itself is explored through scenarios; they note that other, iterative methods can include it at greater computational cost. The Cholesky factor is triangular and is not the symmetric matrix square root, which Higham distinguishes from it. A covariance matrix copied from a publication with few decimal places may fail to factorise even though the original did not, and the rows of the factor follow the order of the covariance matrix, so sampled values must be mapped back to the same parameters in the model.

Sources

• Briggs AH, Weinstein MC, Fenwick EAL, Karnon J, Sculpher MJ, Paltiel AD. Model parameter estimation and uncertainty: a report of the ISPOR-SMDM Modeling Good Research Practices Task Force-6. Value in Health. 2012;15(6):835-842.
• Bujkiewicz S, Achana F, Papanikos T, Riley RD, Abrams KR. NICE DSU Technical Support Document 20: Multivariate meta-analysis of summary data for combining treatment effects on correlated outcomes and evaluating surrogate endpoints. Sheffield: Decision Support Unit, ScHARR; 2019, version 2 updated December 2022. Section 4.1 and Appendix C.1.
• CADTH. Guidelines for the Economic Evaluation of Health Technologies: Canada. 4th ed. Ottawa: CADTH; 2017. Section 12.4 and probabilistic analysis guidance.
• Higham NJ. Cholesky factorization. WIREs Computational Statistics. 2009;1(2):251-254.
• Higham NJ. Computing the nearest correlation matrix: a problem from finance. IMA Journal of Numerical Analysis. 2002;22(3):329-343.
• Latimer N. NICE DSU Technical Support Document 14: Survival analysis for economic evaluations alongside clinical trials, extrapolation with patient-level data. Sheffield: Decision Support Unit, ScHARR; 2011, updated March 2013. Section 3.7.
• LAPACK. DPOTRF reference implementation. Netlib; accessed 9 October 2026. https://www.netlib.org/lapack/double/dpotrf.f
• Neine M, Curran D. An algorithm to generate correlated input-parameters to be used in probabilistic sensitivity analyses. Journal of Market Access and Health Policy. 2021;9(1):1857052.
• NumPy. numpy.linalg.cholesky. NumPy v2.5 manual; accessed 9 October 2026. https://numpy.org/doc/stable/reference/generated/numpy.linalg.cholesky.html
• R Core Team. chol: The Cholesky decomposition. R documentation, base package; accessed 9 October 2026. https://stat.ethz.ch/R-manual/R-devel/library/base/html/chol.html
• Zorginstituut Nederland. Guideline for economic evaluations in healthcare. 2024 version. Diemen: Zorginstituut Nederland; 2024.