<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Implement estimators of large-scale sparse Gaussian densities
by Soumyajit De (email
Step1:  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 <a href="http
Step2:  Next, we use <a href="http
Step3:  <p>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.</p>
<p>Then we define <a href="http
Step4:  Finally, we use the <a href="http
Step5:  To verify the accuracy of the estimate, we compute exact log-determinant of A using Cholesky factorization using <a href="http
Step6:  <h2>Statistics</h2>
We use a smaller sparse-matrix, <a href="http
Step7:  <h2>A motivational example - likelihood of the Ozone dataset</h2>
<p>In <a href="http
Step8:  <h2>Useful components</h2>
<p>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.
<h3>1. <a href="http
Step9:  <h3>2. <a href="http
Step10:  <h4><a href="http
Step11:  <h4><a href="http
Step12:  Apart from iterative solvers, a few more triangular solvers are added.
<h4><a href="http
Step13:  <h4><a href="http
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
from scipy.sparse import eye
from scipy.io import mmread
from matplotlib import pyplot as plt
import os
SHOGUN_DATA_DIR=os.getenv('SHOGUN_DATA_DIR', '../../../data')
matFile=os.path.join(SHOGUN_DATA_DIR, 'logdet/apache2.mtx.gz')
M = mmread(matFile)
rows = M.shape[0]
cols = M.shape[1]
A = M + eye(rows, cols) * 10000.0
plt.title("A")
plt.spy(A, precision = 1e-2, marker = '.', markersize = 0.01)
plt.show()
Explanation: Implement estimators of large-scale sparse Gaussian densities
by Soumyajit De (email: heavensdevil6909@gmail.com, soumyajitde@cse.iitb.ac.in. Github: <a href="https://github.com/lambday">lambday</a>)<br/> Many many thanks to my mentor Heiko Strathmann, Sergey Lisitsyn, Sören Sonnenburg, Viktor Gal
This notebook illustrates large-scale sparse Gaussian density likelihood estimation. It first introduces the reader to the mathematical background and then shows how one can do the estimation with Shogun on a number of real-world data sets.
<h2>Theoretical introduction</h2>
<p><i>Multivariate Gaussian distributions</i>, i.e. some random vector  $\mathbf{x}\in\mathbb{R}^n$ having probability density function
$$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)$$
$\boldsymbol\mu$ being the mean vector and $\boldsymbol\Sigma$ being the covariance matrix, arise in numerous occassions involving large datasets. Computing <i>log-likelihood</i> in these requires computation of the log-determinant of the covariance matrix
$$\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)$$
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$.</p>
<p>The usual technique for computing the log-determinant term in the likelihood expression relies on <i><a href="http://en.wikipedia.org/wiki/Cholesky_factorization">Cholesky factorization</a></i> 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 <i>fill-in</i> 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. <i><a href="http://en.wikipedia.org/wiki/Minimum_degree_algorithm">approximate minimum degree</a></i> (AMD) reordering, these techniques depend largely on the sparsity pattern and therefore not guaranteed to give better result.</p>
<p>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 [<a href="http://link.springer.com/article/10.1007%2Fs11222-012-9368-y">Aune et. al., 2012</a>]. 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
$$\text{trace}(\log(\boldsymbol\Sigma))=\sum_{j=1}^{n}\mathbf{e}^{T}_{j}\log(\boldsymbol\Sigma)\mathbf{e}_{j}$$
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 <i><a href="http://en.wikipedia.org/wiki/Monte_Carlo_method">Monte Carlo estimator</a></i> for the trace -
$$\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}$$
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 <i>probing-vectors</i> that are obtained from <a href="http://en.wikipedia.org/wiki/Graph_coloring">coloring of the adjacency graph</a> represented by the power of the sparse-matrix, $\boldsymbol\Sigma^{p}$, i.e. we can obtain -
$$\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}$$
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.</p>
<p>There are two main issues in this technique. First, computing $\boldsymbol\Sigma^{p}$ is computationally costly, but experiments show that directly applying a <i>d-distance</i> 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 <a href="http://eprints.ma.man.ac.uk/1136/01/covered/MIMS_ep2007_103.pdf">Hale et. al [2008]</a>, i.e. writing $\log(\boldsymbol\Sigma)\mathbf{w}_{j}$ in our desired expression using <i><a href="http://en.wikipedia.org/wiki/Cauchy's_integral_formula">Cauchy's integral formula</a></i> as -
$$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)$$
$K$, $k \in \mathbb{R}$,  $\alpha_{l}$, $\sigma_{l} \in \mathbb{C}$ are coming from <i><a href="http://en.wikipedia.org/wiki/Jacobi_elliptic_functions">Jacobi elliptic functions</a></i>, $\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}$.</p>
<p>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 <i>Krylov subspace iterative solvers</i> like <i><a href="http://en.wikipedia.org/wiki/Conjugate_gradient_method">Conjugate Gradient</a></i> (CG). Since the shifted matrices have complex entries along their diagonal, the appropriate method to choose is <i>Conjugate Orthogonal Conjugate Gradient</i> (COCG) [<a href="http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=106415&tag=1">H.A. van der Vorst et. al., 1990.</a>]. Alternatively, these systems can be solved at once using <i>CG-M</i> [<a href"http://arxiv.org/abs/hep-lat/9612014">Jegerlehner, 1996.</a>] 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.</p>
<p>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.</p>
<h2>An example of using this approach in Shogun</h2>
<p>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, <a href="http://www.cise.ufl.edu/research/sparse/matrices/GHS_psdef/apache2.html">apache2</a>, which is obtained from the <a href="http://www.cise.ufl.edu/research/sparse/matrices/">The University of Florida Sparse Matrix Collection</a>. 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. 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.
End of explanation
from shogun import RealSparseMatrixOperator, LanczosEigenSolver
op = RealSparseMatrixOperator(A.tocsc())
# Lanczos iterative Eigensolver to compute the min/max Eigenvalues which is required to compute the shifts
eigen_solver = LanczosEigenSolver(op)
# we set the iteration limit high to compute the eigenvalues more accurately, default iteration limit is 1000
eigen_solver.set_max_iteration_limit(2000)
# computing the eigenvalues
eigen_solver.compute()
print('Minimum Eigenvalue:', eigen_solver.get_min_eigenvalue())
print('Maximum Eigenvalue:', eigen_solver.get_max_eigenvalue())
Explanation: 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 <a href="http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1SparseMatrixOperator.html">RealSparseMatrixOperator</a> that is suitable for this example. All the solvers work with <a href="http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1LinearOperator.html">LinearOperator</a> type objects. For computing the eigenvalues, we use <a href="http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1LanczosEigenSolver.html">LanczosEigenSolver</a> class. Although computation of the Eigenvalues is done internally within the log-determinant estimator itself (see below), here we explicitely precompute them.
End of explanation
# We can specify the power of the sparse-matrix that is to be used for coloring, default values will apply a
# 2-distance greedy graph coloring algorithm on the sparse-matrix itself. Matrix-power, if specified, is computed in O(lg p)
from shogun import ProbingSampler
trace_sampler = ProbingSampler(op)
# apply the graph coloring algorithm and generate the number of colors, i.e. number of trace samples
trace_sampler.precompute()
print('Number of colors used:', trace_sampler.get_num_samples())
Explanation: Next, we use <a href="http://www.shogun-toolbox.org/doc/en/latest/ProbingSampler_8h_source.html">ProbingSampler</a> class which uses an external library <a href="http://www.cscapes.org/coloringpage/">ColPack</a>. Again, the number of colors used is precomputed for demonstration purpose, although computed internally inside the log-determinant estimator.
End of explanation
from shogun import CGMShiftedFamilySolver, LogRationalApproximationCGM
cgm = CGMShiftedFamilySolver()
# setting the iteration limit (set this to higher value for higher condition number)
cgm.set_iteration_limit(100)
# accuracy determines the number of contour points in the rational approximation (i.e. number of shifts in the systems)
accuracy = 1E-15
# we create a operator-log-function using the sparse matrix operator that uses CG-M to solve the shifted systems
op_func = LogRationalApproximationCGM(op, eigen_solver, cgm, accuracy)
op_func.precompute()
print('Number of shifts:', op_func.get_num_shifts())
Explanation: <p>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.</p>
<p>Then we define <a href="http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CLogRationalApproximationCGM.html">LogRationalApproximationCGM</a> operator function class, which internally uses the Eigensolver to compute the Eigenvalues, uses <a href="http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CJacobiEllipticFunctions.html">JacobiEllipticFunctions</a> 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 then uses CG-M solver (<a href="http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CGMShiftedFamilySolver.html">CGMShiftedFamilySolver</a>) to solve the shifted systems. Precompute is not necessary here too.</p>
End of explanation
import numpy as np
from shogun import LogDetEstimator
# number of log-det samples (use a higher number to get better estimates)
# (this is 5 times number of colors estimate in practice, so usually 1 probing estimate is enough)
num_samples = 5
log_det_estimator = LogDetEstimator(trace_sampler, op_func)
estimates = log_det_estimator.sample(num_samples)
estimated_logdet = np.mean(estimates)
print('Estimated log(det(A)):', estimated_logdet)
Explanation: Finally, we use the <a href="http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1LogDetEstimator.html">LogDetEstimator</a> class to sample the log-determinant of the matrix.
End of explanation
# the following method requires massive amount of memory, for demonstration purpose
# the following code is commented out and direct value obtained from running it once is used
# from shogun import Statistics
# actual_logdet = Statistics.log_det(A)
actual_logdet = 7120357.73878
print('Actual log(det(A)):', actual_logdet)
plt.hist(estimates)
plt.plot([actual_logdet, actual_logdet], [0,len(estimates)], linewidth=3)
plt.show()
Explanation: To verify the accuracy of the estimate, we compute exact log-determinant of A using Cholesky factorization using <a href="http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1Statistics.html#a9931a4ea72310b239efdc05503442525">Statistics::log_det</a> method.
End of explanation
from scipy.sparse import csc_matrix
from scipy.sparse import identity
m = mmread(os.path.join(SHOGUN_DATA_DIR, 'logdet/west0479.mtx'))
# computing a spd with added ridge
B = csc_matrix(m.transpose() * m + identity(m.shape[0]) * 1000.0)
fig = plt.figure(figsize=(12, 4))
ax = fig.add_subplot(1,2,1)
ax.set_title('B')
ax.spy(B, precision = 1e-5, marker = '.', markersize = 2.0)
ax = fig.add_subplot(1,2,2)
ax.set_title('lower Cholesky factor')
dense_matrix = B.todense()
L = np.linalg.cholesky(dense_matrix)
ax.spy(csc_matrix(L), precision = 1e-5, marker = '.', markersize = 2.0)
plt.show()
op = RealSparseMatrixOperator(B)
eigen_solver = LanczosEigenSolver(op)
# computing log-det estimates using probing sampler
probing_sampler = ProbingSampler(op)
cgm.set_iteration_limit(500)
op_func = LogRationalApproximationCGM(op, eigen_solver, cgm, 1E-5)
log_det_estimator = LogDetEstimator(probing_sampler, op_func)
num_probing_estimates = 100
probing_estimates = log_det_estimator.sample(num_probing_estimates)
# computing log-det estimates using Gaussian sampler
from shogun import NormalSampler, Statistics
num_colors = probing_sampler.get_num_samples()
normal_sampler = NormalSampler(op.get_dimension())
log_det_estimator = LogDetEstimator(normal_sampler, op_func)
num_normal_estimates = num_probing_estimates * num_colors
normal_estimates = log_det_estimator.sample(num_normal_estimates)
# average in groups of n_effective_samples
effective_estimates_normal = np.zeros(num_probing_estimates)
for i in range(num_probing_estimates):
    idx = i * num_colors
    effective_estimates_normal[i] = np.mean(normal_estimates[idx:(idx + num_colors)])
actual_logdet = Statistics.log_det(B)
print('Actual log(det(B)):', actual_logdet)
print('Estimated log(det(B)) using probing sampler:', np.mean(probing_estimates))
print('Estimated log(det(B)) using Gaussian sampler:', np.mean(effective_estimates_normal))
print('Variance using probing sampler:', np.var(probing_estimates))
print('Variance using Gaussian sampler:', np.var(effective_estimates_normal))
fig = plt.figure(figsize=(15, 4))
ax = fig.add_subplot(1,3,1)
ax.set_title('Probing sampler')
ax.plot(np.cumsum(probing_estimates)/(np.arange(len(probing_estimates))+1))
ax.plot([0,len(probing_estimates)], [actual_logdet, actual_logdet])
ax.legend(["Probing", "True"])
ax = fig.add_subplot(1,3,2)
ax.set_title('Gaussian sampler')
ax.plot(np.cumsum(effective_estimates_normal)/(np.arange(len(effective_estimates_normal))+1))
ax.plot([0,len(probing_estimates)], [actual_logdet, actual_logdet])
ax.legend(["Gaussian", "True"])
ax = fig.add_subplot(1,3,3)
ax.hist(probing_estimates)
ax.hist(effective_estimates_normal)
ax.plot([actual_logdet, actual_logdet], [0,len(probing_estimates)], linewidth=3)
plt.show()
Explanation: <h2>Statistics</h2>
We use a smaller sparse-matrix, <a href="http://www.cise.ufl.edu/research/sparse/matrices/HB/west0479.html">'west0479'</a> in this section to demonstrate the benefits of using probing vectors over standard Gaussian vectors to sample the trace of matrix-logarithm. In the following we can easily observe the fill-in phenomena described earlier. Again, a ridge has been added to reduce the runtime for demonstration purpose.
End of explanation
from scipy.io import loadmat
def get_Q_y_A(kappa):
    # read the ozone data and create the matrix Q
    ozone = loadmat(os.path.join(SHOGUN_DATA_DIR, 'logdet/ozone_data.mat'))
    GiCG = ozone["GiCG"]
    G = ozone["G"]
    C0 = ozone["C0"]
    kappa = 13.1
    Q = GiCG + 2 * (kappa ** 2) * G + (kappa ** 4) * C0
    # also, added a ridge here
    Q = Q + eye(Q.shape[0], Q.shape[1]) * 10000.0
    
    plt.spy(Q, precision = 1e-5, marker = '.', markersize = 1.0)
    plt.show()
    
    # read y and A
    y = ozone["y_ozone"]
    A = ozone["A"]
    return Q, y, A
    
def log_det(A):
    op = RealSparseMatrixOperator(A)
    eigen_solver = LanczosEigenSolver(op)
    probing_sampler = ProbingSampler(op)
    cgm = CGMShiftedFamilySolver()
    cgm.set_iteration_limit(100)
    op_func = LogRationalApproximationCGM(op, eigen_solver, cgm, 1E-5)
    log_det_estimator = LogDetEstimator(probing_sampler, op_func)
    num_estimates = 1
    return np.mean(log_det_estimator.sample(num_estimates))
def log_likelihood(tau, kappa):
    Q, y, A = get_Q_y_A(kappa)
    n = len(y);
    AtA = A.T.dot(A)
    M = Q + tau * AtA;
    # Computing log-determinants")
    logdet1 = log_det(Q)
    logdet2 = log_det(M)
        
    first = 0.5 * logdet1 + 0.5 * n * np.log(tau) - 0.5 * logdet2
    
    # computing the rest of the likelihood
    second_a = -0.5 * tau * (y.T.dot(y))
    
    second_b = np.array(A.T.dot(y))
    from scipy.sparse.linalg import spsolve
    second_b = spsolve(M, second_b)
    second_b = A.dot(second_b)
    second_b = y.T.dot(second_b)
    second_b = 0.5 * (tau ** 2) * second_b
        
    log_det_part = first
    quadratic_part = second_a + second_b
    const_part = -0.5 * n * np.log(2 * np.pi)
      
    log_marignal_lik = const_part + log_det_part + quadratic_part
    return log_marignal_lik
L = log_likelihood(1.0, 15.0)
print('Log-likelihood estimate:', L)
Explanation: <h2>A motivational example - likelihood of the Ozone dataset</h2>
<p>In <a href="http://arxiv.org/abs/1306.4032">Lyne et. al. (2013)</a>, 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 -
$$y_{i}|\mathbf{x},\kappa,\tau\sim\mathcal{N}(\mathbf{Ax},\tau^{−1}\mathbf{I})$$
$$\mathbf{x}|\kappa\sim\mathcal{N}(\mathbf{0}, \mathbf{Q}(\kappa))$$
$$\kappa\sim\log_{2}\mathcal{N}(0, 100), \tau\sim\log_{2}\mathcal{N}(0, 100)$$
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
$$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}$$
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, <a href="https://github.com/karlnapf/ozone-roulette.git">ozone-roulette</a>, written by Heiko Strathmann, one of the authors of the original paper).
<b>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.</b>
End of explanation
from shogun import RealSparseMatrixOperator, ComplexDenseMatrixOperator
dim = 5
np.random.seed(10)
# create a random valued sparse matrix linear operator
A = csc_matrix(np.random.randn(dim, dim))
op = RealSparseMatrixOperator(A)
# creating a random vector
np.random.seed(1)
b = np.array(np.random.randn(dim))
v = op.apply(b)
print('A.apply(b)=',v)
# create a dense matrix linear operator
B = np.array(np.random.randn(dim, dim)).astype(complex)
op = ComplexDenseMatrixOperator(B)
print('Dimension:', op.get_dimension())
Explanation: <h2>Useful components</h2>
<p>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.
<h3>1. <a href="http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1LinearOperator.html">Linear Operators</a></h3>
All the linear solvers and Eigen solvers work with linear operators. Both real valued and complex valued operators are supported for dense/sparse matrix linear operators.
End of explanation
from scipy.sparse import csc_matrix
from scipy.sparse import identity
from shogun import ConjugateGradientSolver
# creating a random spd matrix
dim = 5
np.random.seed(10)
m = csc_matrix(np.random.randn(dim, dim))
a = m.transpose() * m + csc_matrix(np.identity(dim))
Q = RealSparseMatrixOperator(a)
# creating a random vector
y = np.array(np.random.randn(dim))
# solve the system Qx=y
# the argument is set as True to gather convergence statistics (default is False)
cg = ConjugateGradientSolver(True)
cg.set_iteration_limit(20)
x = cg.solve(Q,y)
print('x:',x)
# verifying the result
print('y:', y)
print('Qx:', Q.apply(x))
residuals = cg.get_residuals()
plt.plot(residuals)
plt.show()
Explanation: <h3>2. <a href="http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1LinearSolver.html">Linear Solvers</a></h3>
<p> 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.</p>
<h4> <a href="http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CConjugateGradientSolver.html">Conjugate Gradient Solver</a></h4>
This solver solves for system $\mathbf{Qx}=\mathbf{y}$, where $\mathbf{Q}$ is real-valued spd linear operator (e.g. dense/sparse matrix operator), and $\mathbf{y}$ is real vector.
End of explanation
from shogun import ComplexSparseMatrixOperator
from shogun import ConjugateOrthogonalCGSolver
# creating a random spd matrix
dim = 5
np.random.seed(10)
m = csc_matrix(np.random.randn(dim, dim))
a = m.transpose() * m + csc_matrix(np.identity(dim))
a = a.astype(complex)
# adding a complex entry along the diagonal
for i in range(0, dim):
    a[i,i] += complex(np.random.randn(), np.random.randn())
Q = ComplexSparseMatrixOperator(a)
z = np.array(np.random.randn(dim))
# solve for the system Qx=z
cocg = ConjugateOrthogonalCGSolver(True)
cocg.set_iteration_limit(20)
x = cocg.solve(Q, z)
print('x:',x)
# verifying the result
print('z:',z)
print('Qx:',np.real(Q.apply(x)))
residuals = cocg.get_residuals()
plt.plot(residuals)
plt.show()
Explanation: <h4><a href="http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1ConjugateOrthogonalCGSolver.html">Conjugate Orthogonal CG Solver</a></h4>
Solves for systems $\mathbf{Qx}=\mathbf{z}$, where $\mathbf{Q}$ is symmetric but non-Hermitian (i.e. having complex entries in its diagonal) and $\mathbf{z}$ is real valued vector.
End of explanation
from shogun import CGMShiftedFamilySolver
cgm = CGMShiftedFamilySolver()
# creating a random spd matrix
dim = 5
np.random.seed(10)
m = csc_matrix(np.random.randn(dim, dim))
a = m.transpose() * m + csc_matrix(np.identity(dim))
Q = RealSparseMatrixOperator(a)
# creating a random vector
v = np.array(np.random.randn(dim))
# number of shifts (will be equal to the number of contour points)
num_shifts = 3;
# generating some random shifts
shifts = []
for i in range(0, num_shifts):
    shifts.append(complex(np.random.randn(), np.random.randn()))
sigma = np.array(shifts)
print('Shifts:', sigma)
# generating some random weights
weights = []
for i in range(0, num_shifts):
    weights.append(complex(np.random.randn(), np.random.randn()))
alpha = np.array(weights)
print('Weights:',alpha)
# solve for the systems
cgm = CGMShiftedFamilySolver(True)
cgm.set_iteration_limit(20)
x = cgm.solve_shifted_weighted(Q, v, sigma, alpha)
print('x:',x)
residuals = cgm.get_residuals()
plt.plot(residuals)
plt.show()
# verifying the result with cocg
x_s = np.array([0+0j] * dim)
for i in range(0, num_shifts):
    a_s = a.astype(complex)
    for j in range(0, dim):
        # moving the complex shift inside the operator
        a_s[j,j] += sigma[i]
    Q_s = ComplexSparseMatrixOperator(a_s)
    # multiplying the result with weight
    x_s += alpha[i] * cocg.solve(Q_s, v)
print('x\':', x_s)
Explanation: <h4><a href="http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CGMShiftedFamilySolver.html">CG-M Shifted Family Solver</a></h4>
Solves for systems with real valued spd matrices with complex shifts. For using it with log-det, an option to specify the weight of each solution is also there. The solve_shifted_weighted method returns $\sum\alpha_{l}\mathbf{x}{l}$ where $\mathbf{x}{l}=(\mathbf{A}+\sigma_{l}\mathbf{I})^{-1}\mathbf{y}$, $\sigma,\alpha\in\mathbb{C}$, $\mathbf{y}\in\mathbb{R}$.
End of explanation
from shogun import DirectSparseLinearSolver
# creating a random spd matrix
dim = 5
np.random.seed(10)
m = csc_matrix(np.random.randn(dim, dim))
a = m.transpose() * m + csc_matrix(np.identity(dim))
Q = RealSparseMatrixOperator(a)
# creating a random vector
y = np.array(np.random.randn(dim))
# solve the system Qx=y
chol = DirectSparseLinearSolver()
x = chol.solve(Q,y)
print('x:',x)
# verifying the result
print('y:', y)
print('Qx:', Q.apply(x))
Explanation: Apart from iterative solvers, a few more triangular solvers are added.
<h4><a href="http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CDirectSparseLinearSolver.html">Direct Sparse Linear Solver</a></h4>
This uses sparse Cholesky to solve for linear systems $\mathbf{Qx}=\mathbf{y}$, where $\mathbf{Q}$ is real-valued spd linear operator (e.g. dense/sparse matrix operator), and $\mathbf{y}$ is real vector.
End of explanation
from shogun import DirectLinearSolverComplex
# creating a random spd matrix
dim = 5
np.random.seed(10)
m = np.array(np.random.randn(dim, dim))
a = m.transpose() * m + csc_matrix(np.identity(dim))
a = a.astype(complex)
# adding a complex entry along the diagonal
for i in range(0, dim):
    a[i,i] += complex(np.random.randn(), np.random.randn())
Q = ComplexDenseMatrixOperator(a)
z = np.array(np.random.randn(dim))
# solve for the system Qx=z
solver = DirectLinearSolverComplex()
x = solver.solve(Q, z)
print('x:',x)
# verifying the result
print('z:',z)
print('Qx:',np.real(Q.apply(x)))
Explanation: <h4><a href="http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CDirectLinearSolverComplex.html">Direct Linear Solver for Complex</a></h4>
This solves linear systems $\mathbf{Qx}=\mathbf{y}$, where $\mathbf{Q}$ is complex-valued dense matrix linear operator, and $\mathbf{y}$ is real vector.
End of explanation
<END_TASK>