Skip to content

Commit

Permalink
fix notebook after removal of CComputationEngine framework
Browse files Browse the repository at this point in the history
  • Loading branch information
@karlnapf
karlnapf committed May 22, 2018
1 parent 1e0b537 commit 567fdc6
Showing 1 changed file with 38 additions and 48 deletions.
86 changes: 38 additions & 48 deletions doc/ipython-notebooks/logdet/logdet.ipynb
View file
Expand Up @@ -18,7 +18,7 @@
"version": "2.7.12"
},
"name": "",
"signature": "sha256:285eedd45763f28c1bb30efd2f2f61c3b5a9dbe1c021a246bc2076c5350f8bc7"
"signature": "sha256:77bb141c1b2a113aa221ec383efb5ac3fad3d7f7cdf683a3eba8dfe2d8b0a9ed"
},
"nbformat": 3,
"nbformat_minor": 0,
Expand Down Expand Up @@ -78,7 +78,7 @@
"<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>\n",
"\n",
"<h2>An example of using this approach in Shogun</h2>\n",
"<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 which is then used with <a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CSerialComputationEngine.html\">SerialComputationEngine</a> 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."
"<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."
]
},
{
Expand All @@ -103,8 +103,7 @@
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": null
"outputs": []
},
{
"cell_type": "markdown",
Expand Down Expand Up @@ -133,8 +132,7 @@
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": null
"outputs": []
},
{
"cell_type": "markdown",
Expand All @@ -158,25 +156,23 @@
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": null
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<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>\n",
"\n",
"<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 submits a computation job to the engine which then uses CG-M solver (<a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CCGMShiftedFamilySolver.html\">CGMShiftedFamilySolver</a>) to solve the shifted systems. Precompute is not necessary here too.</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_1CCGMShiftedFamilySolver.html\">CGMShiftedFamilySolver</a>) to solve the shifted systems. Precompute is not necessary here too.</p>"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from shogun import SerialComputationEngine, CGMShiftedFamilySolver, LogRationalApproximationCGM\n",
"from shogun import 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",
Expand All @@ -185,14 +181,13 @@
"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 = LogRationalApproximationCGM(op, 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
"outputs": []
},
{
"cell_type": "markdown",
Expand All @@ -212,16 +207,15 @@
"# (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",
"log_det_estimator = LogDetEstimator(trace_sampler, op_func)\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
"outputs": []
},
{
"cell_type": "markdown",
Expand Down Expand Up @@ -249,8 +243,7 @@
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": null
"outputs": []
},
{
"cell_type": "markdown",
Expand Down Expand Up @@ -286,8 +279,7 @@
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": null
"outputs": []
},
{
"cell_type": "code",
Expand All @@ -300,8 +292,8 @@
"probing_sampler = ProbingSampler(op)\n",
"cgm.set_iteration_limit(500)\n",
"\n",
"op_func = LogRationalApproximationCGM(op, engine, eigen_solver, cgm, 1E-5)\n",
"log_det_estimator = LogDetEstimator(probing_sampler, op_func, engine)\n",
"op_func = LogRationalApproximationCGM(op, eigen_solver, cgm, 1E-5)\n",
"log_det_estimator = LogDetEstimator(probing_sampler, op_func)\n",
"num_probing_estimates = 100\n",
"probing_estimates = log_det_estimator.sample(num_probing_estimates)\n",
"\n",
Expand All @@ -310,7 +302,7 @@
"\n",
"num_colors = probing_sampler.get_num_samples()\n",
"normal_sampler = NormalSampler(op.get_dimension())\n",
"log_det_estimator = LogDetEstimator(normal_sampler, op_func, engine)\n",
"log_det_estimator = LogDetEstimator(normal_sampler, op_func)\n",
"num_normal_estimates = num_probing_estimates * num_colors\n",
"normal_estimates = log_det_estimator.sample(num_normal_estimates)\n",
"\n",
Expand All @@ -321,16 +313,15 @@
" effective_estimates_normal[i] = np.mean(normal_estimates[idx:(idx + num_colors)])\n",
"\n",
"actual_logdet = Statistics.log_det(B)\n",
"print('Actual log(det(B)):', actual_logdet\n",
"print('Actual log(det(B)):', actual_logdet)\n",
"print('Estimated log(det(B)) using probing sampler:', np.mean(probing_estimates))\n",
"print('Estimated log(det(B)) using Gaussian sampler:', np.mean(effective_estimates_normal))\n",
"print('Variance using probing sampler:', np.var(probing_estimates))\n",
"print('Variance using Gaussian sampler:', np.var(effective_estimates_normal))"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": null
"outputs": []
},
{
"cell_type": "code",
Expand All @@ -357,15 +348,14 @@
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": null
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<h2>A motivational example - likelihood of the Ozone dataset</h2>\n",
"<p>In <a href=\"http://arxiv.org/abs/1306.4032\">Girolami 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 -\n",
"<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 -\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",
Expand Down Expand Up @@ -403,13 +393,12 @@
" \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",
" op_func = LogRationalApproximationCGM(op, eigen_solver, cgm, 1E-5)\n",
" log_det_estimator = LogDetEstimator(probing_sampler, op_func)\n",
" num_estimates = 1\n",
" return np.mean(log_det_estimator.sample(num_estimates))\n",
"\n",
Expand Down Expand Up @@ -448,8 +437,7 @@
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": null
"outputs": []
},
{
"cell_type": "markdown",
Expand Down Expand Up @@ -492,8 +480,7 @@
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": null
"outputs": []
},
{
"cell_type": "markdown",
Expand Down Expand Up @@ -541,8 +528,7 @@
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": null
"outputs": []
},
{
"cell_type": "markdown",
Expand Down Expand Up @@ -589,8 +575,7 @@
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": null
"outputs": []
},
{
"cell_type": "markdown",
Expand Down Expand Up @@ -660,8 +645,7 @@
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": null
"outputs": []
},
{
"cell_type": "markdown",
Expand Down Expand Up @@ -700,8 +684,7 @@
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": null
"outputs": []
},
{
"cell_type": "markdown",
Expand Down Expand Up @@ -742,8 +725,7 @@
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": null
"outputs": []
},
{
"cell_type": "markdown",
Expand All @@ -755,8 +737,16 @@
"<li> Nicholas Hale, Nicholas J. Higham and Lloyd N. Trefethen, <i>Computing $A^{\\alpha}$, $\\log(A)$ and Related Matrix Functions by Contour Integrals</i>, MIMS EPrint: 2007.103</li>\n",
"<li> H.A. van der Vorst, <i>A Petrov-Galerkin Type Method for Solving $\\mathbf{Ax}=\\mathbf{b}$ Where $\\mathbf{A}$ Is Symmetric Complex</i>, IEEE TRANSACTIONS ON MAGNETICS, VOL. 26, NO. 2, MARCH 1990</li>\n",
"<li> B. Jegerlehner, <i>Krylov space solvers for shifted linear systems</i>, HEP-LAT heplat/9612014, 1996</li>\n",
"<li> Mark Girolami, Anne-Marie Lyne, Heiko Strathmann, Daniel Simpson, Yves Atchade, <i>Playing Russian Roulette with Intractable Likelihoods</i>,arXiv:1306.4032 June 2013</li>\n"
"<li> Anne-Marie Lyne, Mark Girolami, Yves Atchade, Heiko Strathmann, Daniel Simpson<i>Playing Russian Roulette with Intractable Likelihoods</i>,arXiv:1306.4032 June 2013</li>\n"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": []
}
],
"metadata": {}
Expand Down

0 comments on commit 567fdc6

Please sign in to comment.