Multivariate Gaussian distributions, i.e. some random vector $\\mathbf{x}\\in\\mathbb{R}^n$ having probability density function\n", - "$$p(\\mathbf{x}|\\boldsymbol\\mu, \\boldsymbol\\Sigma)=(2\\pi)^{-n/2}\\text{det}(\\boldsymbol\\Sigma)^{-1/2} \\exp\\left(-\\frac{1}{2}(\\mathbf{x}-\\boldsymbol\\mu)^{T}\\boldsymbol\\Sigma^{-1}(\\mathbf{x}-\\boldsymbol\\mu)\\right)$$\n", - "$\\boldsymbol\\mu$ being the mean vector and $\\boldsymbol\\Sigma$ being the covariance matrix, arise in numerous occassions involving large datasets. Computing log-likelihood in these requires computation of the log-determinant of the covariance matrix\n", - "$$\\mathcal{L}(\\mathbf{x}|\\boldsymbol\\mu,\\boldsymbol\\Sigma)=-\\frac{n}{2}\\log(2\\pi)-\\frac{1}{2}\\log(\\text{det}(\\boldsymbol\\Sigma))-\\frac{1}{2}(\\mathbf{x}-\\boldsymbol\\mu)^{T}\\boldsymbol\\Sigma^{-1}(\\mathbf{x}-\\boldsymbol\\mu)$$\n", - "The covariance matrix and its inverse are symmetric positive definite (spd) and are often sparse, e.g. due to conditional independence properties of Gaussian Markov Random Fields (GMRF). Therefore they can be stored efficiently even for large dimension $n$.
\n", - "\n", - "The usual technique for computing the log-determinant term in the likelihood expression relies on Cholesky factorization of the matrix, i.e. $\\boldsymbol\\Sigma=\\mathbf{LL}^{T}$, ($\\mathbf{L}$ is the lower triangular Cholesky factor) and then using the diagonal entries of the factor to compute $\\log(\\text{det}(\\boldsymbol\\Sigma))=2\\sum_{i=1}^{n}\\log(\\mathbf{L}_{ii})$. However, for sparse matrices, as covariance matrices usually are, the Cholesky factors often suffer from fill-in phenomena - they turn out to be not so sparse themselves. Therefore, for large dimensions this technique becomes infeasible because of a massive memory requirement for storing all these irrelevant non-diagonal co-efficients of the factor. While ordering techniques have been developed to permute the rows and columns beforehand in order to reduce fill-in, e.g. approximate minimum degree (AMD) reordering, these techniques depend largely on the sparsity pattern and therefore not guaranteed to give better result.
\n", - "\n", - "Recent research shows that using a number of techniques from complex analysis, numerical linear algebra and greedy graph coloring, we can, however, approximate the log-determinant up to an arbitrary precision [Aune et. al., 2012]. The main trick lies within the observation that we can write $\\log(\\text{det}(\\boldsymbol\\Sigma))$ as $\\text{trace}(\\log(\\boldsymbol\\Sigma))$, where $\\log(\\boldsymbol\\Sigma)$ is the matrix-logarithm. Computing the log-determinant then requires extracting the trace of the matrix-logarithm as\n", - "$$\\text{trace}(\\log(\\boldsymbol\\Sigma))=\\sum_{j=1}^{n}\\mathbf{e}^{T}_{j}\\log(\\boldsymbol\\Sigma)\\mathbf{e}_{j}$$\n", - "where each $\\mathbf{e}_{j}$ is a unit basis vector having a 1 in its $j^{\\text{th}}$ position while rest are zeros and we assume that we can compute $\\log(\\boldsymbol\\Sigma)\\mathbf{e}_{j}$ (explained later). For large dimension $n$, this approach is still costly, so one needs to rely on sampling the trace. For example, using stochastic vectors we can obtain a Monte Carlo estimator for the trace -\n", - "$$\\text{trace}(\\log(\\boldsymbol\\Sigma))=\\mathbb{E}_{\\mathbf{v}}(\\mathbf{v}^{T}\\log(\\boldsymbol\\Sigma)\\mathbf{v})\\approx \\sum_{j=1}^{k}\\mathbf{s}^{T}_{j}\\log(\\boldsymbol\\Sigma)\\mathbf{s}_{j}$$\n", - "where the source vectors ($\\mathbf{s}_{j}$) have zero mean and unit variance (e.g. $\\mathbf{s}_{j}\\sim\\mathcal{N}(\\mathbf{0}, \\mathbf{I}), \\forall j\\in[1\\cdots k]$). But since this is a Monte Carlo method, we need many many samples to get sufficiently accurate approximation. However, by a method suggested in Aune et. al., we can reduce the number of samples required drastically by using probing-vectors that are obtained from coloring of the adjacency graph represented by the power of the sparse-matrix, $\\boldsymbol\\Sigma^{p}$, i.e. we can obtain -\n", - "$$\\mathbb{E}_{\\mathbf{v}}(\\mathbf{v}^{T}\\log(\\boldsymbol\\Sigma)\\mathbf{v})\\approx \\sum_{j=1}^{m}\\mathbf{w}^{T}_{j}\\log(\\boldsymbol\\Sigma)\\mathbf{w}_{j}$$\n", - "with $m\\ll n$, where $m$ is the number of colors used in the graph coloring. For a particular color $j$, the probing vector $\\mathbb{w}_{j}$ is obtained by filling with $+1$ or $-1$ uniformly randomly for entries corresponding to nodes of the graph colored with $j$, keeping the rest of the entries as zeros. Since the matrix is sparse, the number of colors used is usually very small compared to the dimension $n$, promising the advantage of this approach.
\n", - "\n", - "There are two main issues in this technique. First, computing $\\boldsymbol\\Sigma^{p}$ is computationally costly, but experiments show that directly applying a d-distance coloring algorithm on the sparse matrix itself also results in a pretty good approximation. Second, computing the exact matrix-logarithm is often infeasible because its is not guaranteed to be sparse. Aune et. al. suggested that we can rely on rational approximation of the matrix-logarithm times vector using an approach described in Hale et. al [2008], i.e. writing $\\log(\\boldsymbol\\Sigma)\\mathbf{w}_{j}$ in our desired expression using Cauchy's integral formula as -\n", - "$$log(\\boldsymbol\\Sigma)\\mathbf{w}_{j}=\\frac{1}{2\\pi i}\\oint_{\\Gamma}log(z)(z\\mathbf{I}-\\boldsymbol\\Sigma)^{-1}\\mathbf{w}_{j}dz\\approx \\frac{-8K(\\lambda_{m}\\lambda_{M})^{\\frac{1}{4}}}{k\\pi N} \\boldsymbol\\Sigma\\Im\\left(-\\sum_{l=1}^{N}\\alpha_{l}(\\boldsymbol\\Sigma-\\sigma_{l}\\mathbf{I})^{-1}\\mathbf{w}_{j}\\right)$$\n", - "$K$, $k \\in \\mathbb{R}$, $\\alpha_{l}$, $\\sigma_{l} \\in \\mathbb{C}$ are coming from Jacobi elliptic functions, $\\lambda_{m}$ and $\\lambda_{M}$ are the minimum/maximum eigenvalues of $\\boldsymbol\\Sigma$ (they have to be real-positive), respectively, $N$ is the number of contour points in the quadrature rule of the above integral and $\\Im(\\mathbf{x})$ represents the imaginary part of $\\mathbf{x}\\in\\mathbb{C}^{n}$.
\n", - "\n", - "The problem then finally boils down to solving the shifted family of linear systems $(\\boldsymbol\\Sigma-\\sigma_{l}\\mathbf{I})\\mathbb{x}_{j}=\\mathbb{w}_{j}$. Since $\\boldsymbol\\Sigma$ is sparse, matrix-vector product is not much costly and therefore these systems can be solved with a low memory-requirement using Krylov subspace iterative solvers like Conjugate Gradient (CG). Since the shifted matrices have complex entries along their diagonal, the appropriate method to choose is Conjugate Orthogonal Conjugate Gradient (COCG) [H.A. van der Vorst et. al., 1990.]. Alternatively, these systems can be solved at once using CG-M [Jegerlehner, 1996.] solver which solves for $(\\mathbf{A}+\\sigma\\mathbf{I})\\mathbf{x}=\\mathbf{b}$ for all values of $\\sigma$ using as many matrix-vector products in the CG-iterations as required to solve for one single shifted system. This algorithm shows reliable convergance behavior for systems with reasonable condition number.
\n", - "\n", - "One interesting property of this approach is that once the graph coloring information and shifts/weights are known, all the computation components - solving linear systems, computing final vector-vector product - are independently computable. Therefore, computation can be speeded up using parallel computation of these. To use this, a computation framework for Shogun is developed and the whole log-det computation works on top of it.
\n", - "\n", - "We demonstrate the usage of this technique to estimate log-determinant of a real-valued spd sparse matrix with dimension $715,176\\times 715,176$ with $4,817,870$ non-zero entries, apache2, which is obtained from the The University of Florida Sparse Matrix Collection. Cholesky factorization with AMD for this sparse-matrix gives rise to factors with $353,843,716$ non-zero entries (from source). We use CG-M solver to solve the shifted systems which is then used with SerialComputationEngine to perform the computation jobs sequentially on a single core, that works even on a normal Desktop machine. Since the original matrix is badly conditioned, here we added a ridge along its diagonal to reduce the condition number so that the CG-M solver converges within reasonable time. Please note that for high condition number, the number of iteration has to be set very high." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%matplotlib inline\n", - "from scipy.sparse import eye\n", - "from scipy.io import mmread\n", - "from matplotlib import pyplot as plt\n", - "import os\n", - "SHOGUN_DATA_DIR=os.getenv('SHOGUN_DATA_DIR', '../../../data')\n", - "\n", - "matFile=os.path.join(SHOGUN_DATA_DIR, 'logdet/apache2.mtx.gz')\n", - "M = mmread(matFile)\n", - "rows = M.shape[0]\n", - "cols = M.shape[1]\n", - "A = M + eye(rows, cols) * 10000.0\n", - "plt.title(\"A\")\n", - "plt.spy(A, precision = 1e-2, marker = '.', markersize = 0.01)\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "First, to keep the notion of Krylov subspace, we view the matrix as a linear operator that applies on a vector, resulting a new vector. We use RealSparseMatrixOperator that is suitable for this example. All the solvers work with LinearOperator type objects. For computing the eigenvalues, we use LanczosEigenSolver class. Although computation of the Eigenvalues is done internally within the log-determinant estimator itself (see below), here we explicitely precompute them." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from shogun import RealSparseMatrixOperator, LanczosEigenSolver\n", - "\n", - "op = RealSparseMatrixOperator(A.tocsc())\n", - "\n", - "# Lanczos iterative Eigensolver to compute the min/max Eigenvalues which is required to compute the shifts\n", - "eigen_solver = LanczosEigenSolver(op)\n", - "# we set the iteration limit high to compute the eigenvalues more accurately, default iteration limit is 1000\n", - "eigen_solver.set_max_iteration_limit(2000)\n", - "\n", - "# computing the eigenvalues\n", - "eigen_solver.compute()\n", - "print('Minimum Eigenvalue:', eigen_solver.get_min_eigenvalue())\n", - "print('Maximum Eigenvalue:', eigen_solver.get_max_eigenvalue())" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Next, we use ProbingSampler class which uses an external library ColPack. Again, the number of colors used is precomputed for demonstration purpose, although computed internally inside the log-determinant estimator." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# We can specify the power of the sparse-matrix that is to be used for coloring, default values will apply a\n", - "# 2-distance greedy graph coloring algorithm on the sparse-matrix itself. Matrix-power, if specified, is computed in O(lg p)\n", - "from shogun import ProbingSampler\n", - "\n", - "trace_sampler = ProbingSampler(op)\n", - "# apply the graph coloring algorithm and generate the number of colors, i.e. number of trace samples\n", - "trace_sampler.precompute()\n", - "print('Number of colors used:', trace_sampler.get_num_samples())" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "
This corresponds to averaging over 13 source vectors rather than one (but has much lower variance as using 13 Gaussian source vectors). A comparison between the convergence behavior of using probing sampler and Gaussian sampler is presented later.
\n", - "\n", - "Then we define LogRationalApproximationCGM operator function class, which internally uses the Eigensolver to compute the Eigenvalues, uses JacobiEllipticFunctions to compute the complex shifts, weights and the constant multiplier in the rational approximation expression, takes the probing vector generated by the trace sampler and submits a computation job to the engine which then uses CG-M solver (CGMShiftedFamilySolver) to solve the shifted systems. Precompute is not necessary here too.
" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from shogun import SerialComputationEngine, CGMShiftedFamilySolver, LogRationalApproximationCGM\n", - "\n", - "engine = SerialComputationEngine()\n", - "cgm = CGMShiftedFamilySolver()\n", - "# setting the iteration limit (set this to higher value for higher condition number)\n", - "cgm.set_iteration_limit(100)\n", - "\n", - "# accuracy determines the number of contour points in the rational approximation (i.e. number of shifts in the systems)\n", - "accuracy = 1E-15\n", - "\n", - "# we create a operator-log-function using the sparse matrix operator that uses CG-M to solve the shifted systems\n", - "op_func = LogRationalApproximationCGM(op, engine, eigen_solver, cgm, accuracy)\n", - "op_func.precompute()\n", - "print('Number of shifts:', op_func.get_num_shifts())" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Finally, we use the LogDetEstimator class to sample the log-determinant of the matrix." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "from shogun import LogDetEstimator\n", - "\n", - "# number of log-det samples (use a higher number to get better estimates)\n", - "# (this is 5 times number of colors estimate in practice, so usually 1 probing estimate is enough)\n", - "num_samples = 5\n", - "\n", - "log_det_estimator = LogDetEstimator(trace_sampler, op_func, engine)\n", - "estimates = log_det_estimator.sample(num_samples)\n", - "\n", - "estimated_logdet = np.mean(estimates)\n", - "print('Estimated log(det(A)):', estimated_logdet)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "To verify the accuracy of the estimate, we compute exact log-determinant of A using Cholesky factorization using Statistics::log_det method." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# the following method requires massive amount of memory, for demonstration purpose\n", - "# the following code is commented out and direct value obtained from running it once is used\n", - "\n", - "# from shogun import Statistics\n", - "# actual_logdet = Statistics.log_det(A)\n", - "\n", - "actual_logdet = 7120357.73878\n", - "print('Actual log(det(A)):', actual_logdet)\n", - "\n", - "plt.hist(estimates)\n", - "plt.plot([actual_logdet, actual_logdet], [0,len(estimates)], linewidth=3)\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In Girolami et. al. (2013), an interesting scenario is discussed where the log-likelihood of a model involving large spatial dataset is considered. The data, collected by a satellite consists of $N=173,405$ ozone measurements around the globe. The data is modelled using three stage hierarchical way -\n", - "$$y_{i}|\\mathbf{x},\\kappa,\\tau\\sim\\mathcal{N}(\\mathbf{Ax},\\tau^{−1}\\mathbf{I})$$\n", - "$$\\mathbf{x}|\\kappa\\sim\\mathcal{N}(\\mathbf{0}, \\mathbf{Q}(\\kappa))$$\n", - "$$\\kappa\\sim\\log_{2}\\mathcal{N}(0, 100), \\tau\\sim\\log_{2}\\mathcal{N}(0, 100)$$\n", - "Where the precision matrix, $\\mathbf{Q}$, of a Matern SPDE model, defined on a fixed traingulation of the globe, is sparse and the parameter $\\kappa$ controls for the range at which correlations in the field are effectively zero (see Girolami et. al. for details). The log-likelihood estiamate of the posterior using this model is\n", - "$$2\\mathcal{L}=2\\log \\pi(\\mathbf{y}|\\kappa,\\tau)=C+\\log(\\text{det}(\\mathbf{Q}(\\kappa)))+N\\log(\\tau)−\\log(\\text{det}(\\mathbf{Q}(\\kappa)+\\tau \\mathbf{A}^{T}\\mathbf{A}))− \\tau\\mathbf{y}^{T}\\mathbf{y}+\\tau^{2}\\mathbf{y}^{T}\\mathbf{A}(\\mathbf{Q}(\\kappa)+\\tau\\mathbf{A}^{T}\\mathbf{A})^{−1}\\mathbf{A}^{T}\\mathbf{y}$$\n", - "In the expression, we have two terms involving log-determinant of large sparse matrices. The rational approximation approach described in the previous section can readily be applicable to estimate the log-likelihood. The following computation shows the usage of Shogun's log-determinant estimator for estimating this likelihood (code has been adapted from an open source library, ozone-roulette, written by Heiko Strathmann, one of the authors of the original paper).\n", - "\n", - "Please note that we again added a ridge along the diagonal for faster execution of this example. Since the original matrix is badly conditioned, one needs to set the iteration limits very high for both the Eigen solver and the linear solver in absense of precondioning." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from scipy.io import loadmat\n", - "\n", - "def get_Q_y_A(kappa):\n", - " # read the ozone data and create the matrix Q\n", - " ozone = loadmat(os.path.join(SHOGUN_DATA_DIR, 'logdet/ozone_data.mat'))\n", - " GiCG = ozone[\"GiCG\"]\n", - " G = ozone[\"G\"]\n", - " C0 = ozone[\"C0\"]\n", - " kappa = 13.1\n", - " Q = GiCG + 2 * (kappa ** 2) * G + (kappa ** 4) * C0\n", - " # also, added a ridge here\n", - " Q = Q + eye(Q.shape[0], Q.shape[1]) * 10000.0\n", - " \n", - " plt.spy(Q, precision = 1e-5, marker = '.', markersize = 1.0)\n", - " plt.show()\n", - " \n", - " # read y and A\n", - " y = ozone[\"y_ozone\"]\n", - " A = ozone[\"A\"]\n", - " return Q, y, A\n", - " \n", - "def log_det(A):\n", - " op = RealSparseMatrixOperator(A)\n", - " engine = SerialComputationEngine()\n", - " eigen_solver = LanczosEigenSolver(op)\n", - " probing_sampler = ProbingSampler(op)\n", - " cgm = CGMShiftedFamilySolver()\n", - " cgm.set_iteration_limit(100)\n", - " op_func = LogRationalApproximationCGM(op, engine, eigen_solver, cgm, 1E-5)\n", - " log_det_estimator = LogDetEstimator(probing_sampler, op_func, engine)\n", - " num_estimates = 1\n", - " return np.mean(log_det_estimator.sample(num_estimates))\n", - "\n", - "def log_likelihood(tau, kappa):\n", - " Q, y, A = get_Q_y_A(kappa)\n", - " n = len(y);\n", - " AtA = A.T.dot(A)\n", - " M = Q + tau * AtA;\n", - "\n", - " # Computing log-determinants\")\n", - " logdet1 = log_det(Q)\n", - " logdet2 = log_det(M)\n", - " \n", - " first = 0.5 * logdet1 + 0.5 * n * np.log(tau) - 0.5 * logdet2\n", - " \n", - " # computing the rest of the likelihood\n", - " second_a = -0.5 * tau * (y.T.dot(y))\n", - " \n", - " second_b = np.array(A.T.dot(y))\n", - " from scipy.sparse.linalg import spsolve\n", - " second_b = spsolve(M, second_b)\n", - " second_b = A.dot(second_b)\n", - " second_b = y.T.dot(second_b)\n", - " second_b = 0.5 * (tau ** 2) * second_b\n", - " \n", - " log_det_part = first\n", - " quadratic_part = second_a + second_b\n", - " const_part = -0.5 * n * np.log(2 * np.pi)\n", - " \n", - " log_marignal_lik = const_part + log_det_part + quadratic_part\n", - "\n", - " return log_marignal_lik\n", - "\n", - "L = log_likelihood(1.0, 15.0)\n", - "print('Log-likelihood estimate:', L)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "
As a part of the implementation of log-determinant estimator, a number of classes have been developed, which may come useful for several other occassions as well." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "
Conjugate Gradient based iterative solvers, that construct the Krylov subspace in their iteration by computing matrix-vector products are most useful for solving sparse linear systems. Here is an overview of CG based solvers that are currently available in Shogun.
\n", - "Multivariate Gaussian distributions, i.e. some random vector $\\mathbf{x}\\in\\mathbb{R}^n$ having probability density function\n", + "$$p(\\mathbf{x}|\\boldsymbol\\mu, \\boldsymbol\\Sigma)=(2\\pi)^{-n/2}\\text{det}(\\boldsymbol\\Sigma)^{-1/2} \\exp\\left(-\\frac{1}{2}(\\mathbf{x}-\\boldsymbol\\mu)^{T}\\boldsymbol\\Sigma^{-1}(\\mathbf{x}-\\boldsymbol\\mu)\\right)$$\n", + "$\\boldsymbol\\mu$ being the mean vector and $\\boldsymbol\\Sigma$ being the covariance matrix, arise in numerous occassions involving large datasets. Computing log-likelihood in these requires computation of the log-determinant of the covariance matrix\n", + "$$\\mathcal{L}(\\mathbf{x}|\\boldsymbol\\mu,\\boldsymbol\\Sigma)=-\\frac{n}{2}\\log(2\\pi)-\\frac{1}{2}\\log(\\text{det}(\\boldsymbol\\Sigma))-\\frac{1}{2}(\\mathbf{x}-\\boldsymbol\\mu)^{T}\\boldsymbol\\Sigma^{-1}(\\mathbf{x}-\\boldsymbol\\mu)$$\n", + "The covariance matrix and its inverse are symmetric positive definite (spd) and are often sparse, e.g. due to conditional independence properties of Gaussian Markov Random Fields (GMRF). Therefore they can be stored efficiently even for large dimension $n$.
\n", + "\n", + "The usual technique for computing the log-determinant term in the likelihood expression relies on Cholesky factorization of the matrix, i.e. $\\boldsymbol\\Sigma=\\mathbf{LL}^{T}$, ($\\mathbf{L}$ is the lower triangular Cholesky factor) and then using the diagonal entries of the factor to compute $\\log(\\text{det}(\\boldsymbol\\Sigma))=2\\sum_{i=1}^{n}\\log(\\mathbf{L}_{ii})$. However, for sparse matrices, as covariance matrices usually are, the Cholesky factors often suffer from fill-in phenomena - they turn out to be not so sparse themselves. Therefore, for large dimensions this technique becomes infeasible because of a massive memory requirement for storing all these irrelevant non-diagonal co-efficients of the factor. While ordering techniques have been developed to permute the rows and columns beforehand in order to reduce fill-in, e.g. approximate minimum degree (AMD) reordering, these techniques depend largely on the sparsity pattern and therefore not guaranteed to give better result.
\n", + "\n", + "Recent research shows that using a number of techniques from complex analysis, numerical linear algebra and greedy graph coloring, we can, however, approximate the log-determinant up to an arbitrary precision [Aune et. al., 2012]. The main trick lies within the observation that we can write $\\log(\\text{det}(\\boldsymbol\\Sigma))$ as $\\text{trace}(\\log(\\boldsymbol\\Sigma))$, where $\\log(\\boldsymbol\\Sigma)$ is the matrix-logarithm. Computing the log-determinant then requires extracting the trace of the matrix-logarithm as\n", + "$$\\text{trace}(\\log(\\boldsymbol\\Sigma))=\\sum_{j=1}^{n}\\mathbf{e}^{T}_{j}\\log(\\boldsymbol\\Sigma)\\mathbf{e}_{j}$$\n", + "where each $\\mathbf{e}_{j}$ is a unit basis vector having a 1 in its $j^{\\text{th}}$ position while rest are zeros and we assume that we can compute $\\log(\\boldsymbol\\Sigma)\\mathbf{e}_{j}$ (explained later). For large dimension $n$, this approach is still costly, so one needs to rely on sampling the trace. For example, using stochastic vectors we can obtain a Monte Carlo estimator for the trace -\n", + "$$\\text{trace}(\\log(\\boldsymbol\\Sigma))=\\mathbb{E}_{\\mathbf{v}}(\\mathbf{v}^{T}\\log(\\boldsymbol\\Sigma)\\mathbf{v})\\approx \\sum_{j=1}^{k}\\mathbf{s}^{T}_{j}\\log(\\boldsymbol\\Sigma)\\mathbf{s}_{j}$$\n", + "where the source vectors ($\\mathbf{s}_{j}$) have zero mean and unit variance (e.g. $\\mathbf{s}_{j}\\sim\\mathcal{N}(\\mathbf{0}, \\mathbf{I}), \\forall j\\in[1\\cdots k]$). But since this is a Monte Carlo method, we need many many samples to get sufficiently accurate approximation. However, by a method suggested in Aune et. al., we can reduce the number of samples required drastically by using probing-vectors that are obtained from coloring of the adjacency graph represented by the power of the sparse-matrix, $\\boldsymbol\\Sigma^{p}$, i.e. we can obtain -\n", + "$$\\mathbb{E}_{\\mathbf{v}}(\\mathbf{v}^{T}\\log(\\boldsymbol\\Sigma)\\mathbf{v})\\approx \\sum_{j=1}^{m}\\mathbf{w}^{T}_{j}\\log(\\boldsymbol\\Sigma)\\mathbf{w}_{j}$$\n", + "with $m\\ll n$, where $m$ is the number of colors used in the graph coloring. For a particular color $j$, the probing vector $\\mathbb{w}_{j}$ is obtained by filling with $+1$ or $-1$ uniformly randomly for entries corresponding to nodes of the graph colored with $j$, keeping the rest of the entries as zeros. Since the matrix is sparse, the number of colors used is usually very small compared to the dimension $n$, promising the advantage of this approach.
\n", + "\n", + "There are two main issues in this technique. First, computing $\\boldsymbol\\Sigma^{p}$ is computationally costly, but experiments show that directly applying a d-distance coloring algorithm on the sparse matrix itself also results in a pretty good approximation. Second, computing the exact matrix-logarithm is often infeasible because its is not guaranteed to be sparse. Aune et. al. suggested that we can rely on rational approximation of the matrix-logarithm times vector using an approach described in Hale et. al [2008], i.e. writing $\\log(\\boldsymbol\\Sigma)\\mathbf{w}_{j}$ in our desired expression using Cauchy's integral formula as -\n", + "$$log(\\boldsymbol\\Sigma)\\mathbf{w}_{j}=\\frac{1}{2\\pi i}\\oint_{\\Gamma}log(z)(z\\mathbf{I}-\\boldsymbol\\Sigma)^{-1}\\mathbf{w}_{j}dz\\approx \\frac{-8K(\\lambda_{m}\\lambda_{M})^{\\frac{1}{4}}}{k\\pi N} \\boldsymbol\\Sigma\\Im\\left(-\\sum_{l=1}^{N}\\alpha_{l}(\\boldsymbol\\Sigma-\\sigma_{l}\\mathbf{I})^{-1}\\mathbf{w}_{j}\\right)$$\n", + "$K$, $k \\in \\mathbb{R}$, $\\alpha_{l}$, $\\sigma_{l} \\in \\mathbb{C}$ are coming from Jacobi elliptic functions, $\\lambda_{m}$ and $\\lambda_{M}$ are the minimum/maximum eigenvalues of $\\boldsymbol\\Sigma$ (they have to be real-positive), respectively, $N$ is the number of contour points in the quadrature rule of the above integral and $\\Im(\\mathbf{x})$ represents the imaginary part of $\\mathbf{x}\\in\\mathbb{C}^{n}$.
\n", + "\n", + "The problem then finally boils down to solving the shifted family of linear systems $(\\boldsymbol\\Sigma-\\sigma_{l}\\mathbf{I})\\mathbb{x}_{j}=\\mathbb{w}_{j}$. Since $\\boldsymbol\\Sigma$ is sparse, matrix-vector product is not much costly and therefore these systems can be solved with a low memory-requirement using Krylov subspace iterative solvers like Conjugate Gradient (CG). Since the shifted matrices have complex entries along their diagonal, the appropriate method to choose is Conjugate Orthogonal Conjugate Gradient (COCG) [H.A. van der Vorst et. al., 1990.]. Alternatively, these systems can be solved at once using CG-M [Jegerlehner, 1996.] solver which solves for $(\\mathbf{A}+\\sigma\\mathbf{I})\\mathbf{x}=\\mathbf{b}$ for all values of $\\sigma$ using as many matrix-vector products in the CG-iterations as required to solve for one single shifted system. This algorithm shows reliable convergance behavior for systems with reasonable condition number.
\n", + "\n", + "One interesting property of this approach is that once the graph coloring information and shifts/weights are known, all the computation components - solving linear systems, computing final vector-vector product - are independently computable. Therefore, computation can be speeded up using parallel computation of these. To use this, a computation framework for Shogun is developed and the whole log-det computation works on top of it.
\n", + "\n", + "We demonstrate the usage of this technique to estimate log-determinant of a real-valued spd sparse matrix with dimension $715,176\\times 715,176$ with $4,817,870$ non-zero entries, apache2, which is obtained from the The University of Florida Sparse Matrix Collection. Cholesky factorization with AMD for this sparse-matrix gives rise to factors with $353,843,716$ non-zero entries (from source). We use CG-M solver to solve the shifted systems which is then used with SerialComputationEngine to perform the computation jobs sequentially on a single core, that works even on a normal Desktop machine. Since the original matrix is badly conditioned, here we added a ridge along its diagonal to reduce the condition number so that the CG-M solver converges within reasonable time. Please note that for high condition number, the number of iteration has to be set very high." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "%matplotlib inline\n", + "from scipy.sparse import eye\n", + "from scipy.io import mmread\n", + "from matplotlib import pyplot as plt\n", + "import os\n", + "SHOGUN_DATA_DIR=os.getenv('SHOGUN_DATA_DIR', '../../../data')\n", + "\n", + "matFile=os.path.join(SHOGUN_DATA_DIR, 'logdet/apache2.mtx.gz')\n", + "M = mmread(matFile)\n", + "rows = M.shape[0]\n", + "cols = M.shape[1]\n", + "A = M + eye(rows, cols) * 10000.0\n", + "plt.title(\"A\")\n", + "plt.spy(A, precision = 1e-2, marker = '.', markersize = 0.01)\n", + "plt.show()" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, to keep the notion of Krylov subspace, we view the matrix as a linear operator that applies on a vector, resulting a new vector. We use RealSparseMatrixOperator that is suitable for this example. All the solvers work with LinearOperator type objects. For computing the eigenvalues, we use LanczosEigenSolver class. Although computation of the Eigenvalues is done internally within the log-determinant estimator itself (see below), here we explicitely precompute them." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from shogun import RealSparseMatrixOperator, LanczosEigenSolver\n", + "\n", + "op = RealSparseMatrixOperator(A.tocsc())\n", + "\n", + "# Lanczos iterative Eigensolver to compute the min/max Eigenvalues which is required to compute the shifts\n", + "eigen_solver = LanczosEigenSolver(op)\n", + "# we set the iteration limit high to compute the eigenvalues more accurately, default iteration limit is 1000\n", + "eigen_solver.set_max_iteration_limit(2000)\n", + "\n", + "# computing the eigenvalues\n", + "eigen_solver.compute()\n", + "print('Minimum Eigenvalue:', eigen_solver.get_min_eigenvalue())\n", + "print('Maximum Eigenvalue:', eigen_solver.get_max_eigenvalue())" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we use ProbingSampler class which uses an external library ColPack. Again, the number of colors used is precomputed for demonstration purpose, although computed internally inside the log-determinant estimator." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# We can specify the power of the sparse-matrix that is to be used for coloring, default values will apply a\n", + "# 2-distance greedy graph coloring algorithm on the sparse-matrix itself. Matrix-power, if specified, is computed in O(lg p)\n", + "from shogun import ProbingSampler\n", + "\n", + "trace_sampler = ProbingSampler(op)\n", + "# apply the graph coloring algorithm and generate the number of colors, i.e. number of trace samples\n", + "trace_sampler.precompute()\n", + "print('Number of colors used:', trace_sampler.get_num_samples())" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
This corresponds to averaging over 13 source vectors rather than one (but has much lower variance as using 13 Gaussian source vectors). A comparison between the convergence behavior of using probing sampler and Gaussian sampler is presented later.
\n", + "\n", + "Then we define LogRationalApproximationCGM operator function class, which internally uses the Eigensolver to compute the Eigenvalues, uses JacobiEllipticFunctions to compute the complex shifts, weights and the constant multiplier in the rational approximation expression, takes the probing vector generated by the trace sampler and submits a computation job to the engine which then uses CG-M solver (CGMShiftedFamilySolver) to solve the shifted systems. Precompute is not necessary here too.
" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from shogun import SerialComputationEngine, CGMShiftedFamilySolver, LogRationalApproximationCGM\n", + "\n", + "engine = SerialComputationEngine()\n", + "cgm = CGMShiftedFamilySolver()\n", + "# setting the iteration limit (set this to higher value for higher condition number)\n", + "cgm.set_iteration_limit(100)\n", + "\n", + "# accuracy determines the number of contour points in the rational approximation (i.e. number of shifts in the systems)\n", + "accuracy = 1E-15\n", + "\n", + "# we create a operator-log-function using the sparse matrix operator that uses CG-M to solve the shifted systems\n", + "op_func = LogRationalApproximationCGM(op, engine, eigen_solver, cgm, accuracy)\n", + "op_func.precompute()\n", + "print('Number of shifts:', op_func.get_num_shifts())" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, we use the LogDetEstimator class to sample the log-determinant of the matrix." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "import numpy as np\n", + "from shogun import LogDetEstimator\n", + "\n", + "# number of log-det samples (use a higher number to get better estimates)\n", + "# (this is 5 times number of colors estimate in practice, so usually 1 probing estimate is enough)\n", + "num_samples = 5\n", + "\n", + "log_det_estimator = LogDetEstimator(trace_sampler, op_func, engine)\n", + "estimates = log_det_estimator.sample(num_samples)\n", + "\n", + "estimated_logdet = np.mean(estimates)\n", + "print('Estimated log(det(A)):', estimated_logdet)" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To verify the accuracy of the estimate, we compute exact log-determinant of A using Cholesky factorization using Statistics::log_det method." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# the following method requires massive amount of memory, for demonstration purpose\n", + "# the following code is commented out and direct value obtained from running it once is used\n", + "\n", + "# from shogun import Statistics\n", + "# actual_logdet = Statistics.log_det(A)\n", + "\n", + "actual_logdet = 7120357.73878\n", + "print('Actual log(det(A)):', actual_logdet)\n", + "\n", + "plt.hist(estimates)\n", + "plt.plot([actual_logdet, actual_logdet], [0,len(estimates)], linewidth=3)\n", + "plt.show()" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In Girolami et. al. (2013), an interesting scenario is discussed where the log-likelihood of a model involving large spatial dataset is considered. The data, collected by a satellite consists of $N=173,405$ ozone measurements around the globe. The data is modelled using three stage hierarchical way -\n", + "$$y_{i}|\\mathbf{x},\\kappa,\\tau\\sim\\mathcal{N}(\\mathbf{Ax},\\tau^{\u22121}\\mathbf{I})$$\n", + "$$\\mathbf{x}|\\kappa\\sim\\mathcal{N}(\\mathbf{0}, \\mathbf{Q}(\\kappa))$$\n", + "$$\\kappa\\sim\\log_{2}\\mathcal{N}(0, 100), \\tau\\sim\\log_{2}\\mathcal{N}(0, 100)$$\n", + "Where the precision matrix, $\\mathbf{Q}$, of a Matern SPDE model, defined on a fixed traingulation of the globe, is sparse and the parameter $\\kappa$ controls for the range at which correlations in the field are effectively zero (see Girolami et. al. for details). The log-likelihood estiamate of the posterior using this model is\n", + "$$2\\mathcal{L}=2\\log \\pi(\\mathbf{y}|\\kappa,\\tau)=C+\\log(\\text{det}(\\mathbf{Q}(\\kappa)))+N\\log(\\tau)\u2212\\log(\\text{det}(\\mathbf{Q}(\\kappa)+\\tau \\mathbf{A}^{T}\\mathbf{A}))\u2212 \\tau\\mathbf{y}^{T}\\mathbf{y}+\\tau^{2}\\mathbf{y}^{T}\\mathbf{A}(\\mathbf{Q}(\\kappa)+\\tau\\mathbf{A}^{T}\\mathbf{A})^{\u22121}\\mathbf{A}^{T}\\mathbf{y}$$\n", + "In the expression, we have two terms involving log-determinant of large sparse matrices. The rational approximation approach described in the previous section can readily be applicable to estimate the log-likelihood. The following computation shows the usage of Shogun's log-determinant estimator for estimating this likelihood (code has been adapted from an open source library, ozone-roulette, written by Heiko Strathmann, one of the authors of the original paper).\n", + "\n", + "Please note that we again added a ridge along the diagonal for faster execution of this example. Since the original matrix is badly conditioned, one needs to set the iteration limits very high for both the Eigen solver and the linear solver in absense of precondioning." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from scipy.io import loadmat\n", + "\n", + "def get_Q_y_A(kappa):\n", + " # read the ozone data and create the matrix Q\n", + " ozone = loadmat(os.path.join(SHOGUN_DATA_DIR, 'logdet/ozone_data.mat'))\n", + " GiCG = ozone[\"GiCG\"]\n", + " G = ozone[\"G\"]\n", + " C0 = ozone[\"C0\"]\n", + " kappa = 13.1\n", + " Q = GiCG + 2 * (kappa ** 2) * G + (kappa ** 4) * C0\n", + " # also, added a ridge here\n", + " Q = Q + eye(Q.shape[0], Q.shape[1]) * 10000.0\n", + " \n", + " plt.spy(Q, precision = 1e-5, marker = '.', markersize = 1.0)\n", + " plt.show()\n", + " \n", + " # read y and A\n", + " y = ozone[\"y_ozone\"]\n", + " A = ozone[\"A\"]\n", + " return Q, y, A\n", + " \n", + "def log_det(A):\n", + " op = RealSparseMatrixOperator(A)\n", + " engine = SerialComputationEngine()\n", + " eigen_solver = LanczosEigenSolver(op)\n", + " probing_sampler = ProbingSampler(op)\n", + " cgm = CGMShiftedFamilySolver()\n", + " cgm.set_iteration_limit(100)\n", + " op_func = LogRationalApproximationCGM(op, engine, eigen_solver, cgm, 1E-5)\n", + " log_det_estimator = LogDetEstimator(probing_sampler, op_func, engine)\n", + " num_estimates = 1\n", + " return np.mean(log_det_estimator.sample(num_estimates))\n", + "\n", + "def log_likelihood(tau, kappa):\n", + " Q, y, A = get_Q_y_A(kappa)\n", + " n = len(y);\n", + " AtA = A.T.dot(A)\n", + " M = Q + tau * AtA;\n", + "\n", + " # Computing log-determinants\")\n", + " logdet1 = log_det(Q)\n", + " logdet2 = log_det(M)\n", + " \n", + " first = 0.5 * logdet1 + 0.5 * n * np.log(tau) - 0.5 * logdet2\n", + " \n", + " # computing the rest of the likelihood\n", + " second_a = -0.5 * tau * (y.T.dot(y))\n", + " \n", + " second_b = np.array(A.T.dot(y))\n", + " from scipy.sparse.linalg import spsolve\n", + " second_b = spsolve(M, second_b)\n", + " second_b = A.dot(second_b)\n", + " second_b = y.T.dot(second_b)\n", + " second_b = 0.5 * (tau ** 2) * second_b\n", + " \n", + " log_det_part = first\n", + " quadratic_part = second_a + second_b\n", + " const_part = -0.5 * n * np.log(2 * np.pi)\n", + " \n", + " log_marignal_lik = const_part + log_det_part + quadratic_part\n", + "\n", + " return log_marignal_lik\n", + "\n", + "L = log_likelihood(1.0, 15.0)\n", + "print('Log-likelihood estimate:', L)" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
As a part of the implementation of log-determinant estimator, a number of classes have been developed, which may come useful for several other occassions as well." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
Conjugate Gradient based iterative solvers, that construct the Krylov subspace in their iteration by computing matrix-vector products are most useful for solving sparse linear systems. Here is an overview of CG based solvers that are currently available in Shogun.
\n", + "