We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
Exact GP regression Lapack Error syev
Training a simple GP regression model produces following runtime error randomly at some seeds: 
`RuntimeError: Lapack Error syev : 3 off-diagonal elements didn't converge to zero at /pytorch/aten/src/TH/generic/THTensorLapack.c:388`.
I have been able to reproduce the error with the following code:

```
import math
import torch
import gpytorch
from matplotlib import pyplot as plt

from torch import optim
from gpytorch.kernels import RBFKernel
from gpytorch.means import ConstantMean
from gpytorch.likelihoods import GaussianLikelihood
from gpytorch.random_variables import GaussianRandomVariable

torch.manual_seed(1)
train_x = torch.FloatTensor([[ 6,  6],
       [ 7,  3],
       [ 7,  8],
       [ 8,  6],
       [10,  0]])

train_y = torch.FloatTensor([2.38922077, 0.35414244, 2.24837906, 1.40895557, 0.68164634])

class ExactGPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        super(ExactGPModel, self).__init__(train_x, train_y, likelihood)
        self.mean_module = ConstantMean(constant_bounds=(-10, 10))
        self.covar_module = RBFKernel(log_lengthscale_bounds=(-5, 5))
    
    def forward(self, x):
        mean_x = self.mean_module(x)
        covar_x = self.covar_module(x)
        return GaussianRandomVariable(mean_x, covar_x)

likelihood = GaussianLikelihood(log_noise_bounds=(-5, 5))
model = ExactGPModel(train_x.data, train_y.data, likelihood)

model.train()
likelihood.train()

optimizer = torch.optim.Adam([
    {'params': model.parameters()},  
], lr=0.1)

mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)

training_iter = 50
for i in range(training_iter):
    optimizer.zero_grad()
    output = model(train_x)
    loss = -mll(output, train_y)
    loss.backward()
    print('Iter %d/%d - Loss: %.3f   log_lengthscale: %.3f   log_noise: %.3f' % (
        i + 1, training_iter, loss.data[0],
        model.covar_module.log_lengthscale.data[0, 0],
        model.likelihood.log_noise.data[0]
    ))
    optimizer.step()


```
I am using the latest version of Gpytorch and Pytorch (0.4.0). The above code is identical to the simple gp regression example provided with this repo with changes in train_x and train_y only. Also, the code works fine for seed = 0. 


--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `gpytorch/utils/linear_cg.py`
Content:
```
1 from __future__ import absolute_import
2 from __future__ import division
3 from __future__ import print_function
4 from __future__ import unicode_literals
5 
6 import torch
7 from .. import settings
8 
9 
10 def _default_preconditioner(x):
11     return x.clone()
12 
13 
14 def linear_cg(
15     matmul_closure,
16     rhs,
17     n_tridiag=0,
18     tolerance=1e-6,
19     eps=1e-20,
20     max_iter=None,
21     max_tridiag_iter=None,
22     initial_guess=None,
23     preconditioner=None,
24 ):
25     """
26     Implements the linear conjugate gradients method for (approximately) solving systems of the form
27 
28         lhs result = rhs
29 
30     for positive definite and symmetric matrices.
31 
32     Args:
33       - matmul_closure - a function which performs a left matrix multiplication with lhs_mat
34       - rhs - the right-hand side of the equation
35       - n_tridiag - returns a tridiagonalization of the first n_tridiag columns of rhs
36       - tolerance - stop the solve when the max residual is less than this
37       - eps - noise to add to prevent division by zero
38       - max_iter - the maximum number of CG iterations
39       - max_tridiag_iter - the maximum size of the tridiagonalization matrix
40       - initial_guess - an initial guess at the solution `result`
41       - precondition_closure - a functions which left-preconditions a supplied vector
42 
43     Returns:
44       result - a solution to the system (if n_tridiag is 0)
45       result, tridiags - a solution to the system, and corresponding tridiagonal matrices (if n_tridiag > 0)
46     """
47     # Unsqueeze, if necesasry
48     is_vector = rhs.ndimension() == 1
49     if is_vector:
50         rhs = rhs.unsqueeze(-1)
51 
52     # Some default arguments
53     if max_iter is None:
54         max_iter = settings.max_cg_iterations.value()
55     if max_tridiag_iter is None:
56         max_tridiag_iter = settings.max_lanczos_quadrature_iterations.value()
57     if initial_guess is None:
58         initial_guess = rhs.new(rhs.size()).zero_()
59     if preconditioner is None:
60         preconditioner = _default_preconditioner
61 
62     # If we are running m CG iterations, we obviously can't get more than m Lanczos coefficients
63     if max_tridiag_iter > max_iter:
64         raise RuntimeError("Getting a tridiagonalization larger than the number of CG iterations run is not possible!")
65 
66     # Check matmul_closure object
67     if torch.is_tensor(matmul_closure):
68         matmul_closure = matmul_closure.matmul
69     elif not callable(matmul_closure):
70         raise RuntimeError("matmul_closure must be a tensor, or a callable object!")
71 
72     # Get some constants
73     n_rows = rhs.size(-2)
74     n_iter = min(max_iter, n_rows)
75     n_tridiag_iter = min(max_tridiag_iter, n_rows)
76 
77     # result <- x_{0}
78     result = initial_guess
79 
80     # residual: residual_{0} = b_vec - lhs x_{0}
81     residual = rhs - matmul_closure(result)
82 
83     # Check for NaNs
84     if not torch.equal(residual, residual):
85         raise RuntimeError("NaNs encounterd when trying to perform matrix-vector multiplication")
86 
87     # Sometime we're lucky and the preconditioner solves the system right away
88     residual_norm = residual.norm(2, dim=-2)
89     if (residual_norm < tolerance).all() and not n_tridiag:
90         n_iter = 0  # Skip the iteration!
91 
92     # Otherwise, let's define precond_residual and curr_conjugate_vec
93     else:
94         # precon_residual{0} = M^-1 residual_{0}
95         precond_residual = preconditioner(residual)
96         curr_conjugate_vec = precond_residual
97         residual_inner_prod = precond_residual.mul(residual).sum(-2, keepdim=True)
98 
99         # Define storage matrices
100         mul_storage = residual.new(residual.size())
101         alpha = residual.new(rhs.size(0), 1, rhs.size(-1)) if rhs.ndimension() == 3 else residual.new(1, rhs.size(-1))
102         beta = alpha.new(alpha.size())
103 
104     # Define tridiagonal matrices, if applicable
105     if n_tridiag:
106         if rhs.ndimension() == 3:
107             t_mat = residual.new(n_tridiag_iter, n_tridiag_iter, rhs.size(0), n_tridiag).zero_()
108             alpha_reciprocal = alpha.new(rhs.size(0), n_tridiag)
109         else:
110             t_mat = residual.new(n_tridiag_iter, n_tridiag_iter, n_tridiag).zero_()
111             alpha_reciprocal = alpha.new(n_tridiag)
112 
113         prev_alpha_reciprocal = alpha.new(alpha_reciprocal.size())
114         prev_beta = alpha.new(alpha_reciprocal.size())
115 
116     # Start the iteration
117     for k in range(n_iter):
118         # Get next alpha
119         # alpha_{k} = (residual_{k-1}^T precon_residual{k-1}) / (p_vec_{k-1}^T mat p_vec_{k-1})
120         mvms = matmul_closure(curr_conjugate_vec)
121         torch.mul(curr_conjugate_vec, mvms, out=mul_storage)
122         torch.sum(mul_storage, -2, keepdim=True, out=alpha)
123         alpha.add_(eps)
124         torch.div(residual_inner_prod, alpha, out=alpha)
125 
126         # Update result
127         # result_{k} = result_{k-1} + alpha_{k} p_vec_{k-1}
128         torch.addcmul(result, alpha, curr_conjugate_vec, out=result)
129 
130         # Update residual
131         # residual_{k} = residual_{k-1} - alpha_{k} mat p_vec_{k-1}
132         torch.addcmul(residual, -1, alpha, mvms, out=residual)
133 
134         # If residual are sufficiently small, then exit loop
135         # Alternatively, exit if this is our last iteration
136         torch.norm(residual, 2, dim=-2, out=residual_norm)
137         if (residual_norm < tolerance).all() and not n_tridiag:
138             break
139 
140         # Update precond_residual
141         # precon_residual{k} = M^-1 residual_{k}
142         precond_residual = preconditioner(residual)
143 
144         # beta_{k} = (precon_residual{k}^T r_vec_{k}) / (precon_residual{k-1}^T r_vec_{k-1})
145         residual_inner_prod.add_(eps)
146         torch.reciprocal(residual_inner_prod, out=beta)
147         torch.mul(residual, precond_residual, out=mul_storage)
148         torch.sum(mul_storage, -2, keepdim=True, out=residual_inner_prod)
149         beta.mul_(residual_inner_prod)
150 
151         # Update curr_conjugate_vec
152         # curr_conjugate_vec_{k} = precon_residual{k} + beta_{k} curr_conjugate_vec_{k-1}
153         curr_conjugate_vec.mul_(beta).add_(precond_residual)
154 
155         # Update tridiagonal matrices, if applicable
156         if n_tridiag and k < n_tridiag_iter:
157             alpha_tridiag = alpha.squeeze_(-2).narrow(-1, 0, n_tridiag)
158             beta_tridiag = beta.squeeze_(-2).narrow(-1, 0, n_tridiag)
159             torch.reciprocal(alpha_tridiag, out=alpha_reciprocal)
160 
161             if k == 0:
162                 t_mat[k, k].copy_(alpha_reciprocal)
163             else:
164                 torch.addcmul(alpha_reciprocal, prev_beta, prev_alpha_reciprocal, out=t_mat[k, k])
165                 torch.mul(prev_beta.sqrt_(), prev_alpha_reciprocal, out=t_mat[k, k - 1])
166                 t_mat[k - 1, k].copy_(t_mat[k, k - 1])
167 
168             prev_alpha_reciprocal.copy_(alpha_reciprocal)
169             prev_beta.copy_(beta_tridiag)
170 
171     if is_vector:
172         result = result.squeeze(-1)
173 
174     if n_tridiag:
175         if rhs.ndimension() == 3:
176             return result, t_mat.permute(3, 2, 0, 1).contiguous()
177         else:
178             return result, t_mat.permute(2, 0, 1).contiguous()
179     else:
180         return result
181
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/gpytorch/utils/linear_cg.py b/gpytorch/utils/linear_cg.py
--- a/gpytorch/utils/linear_cg.py
+++ b/gpytorch/utils/linear_cg.py
@@ -113,6 +113,8 @@
         prev_alpha_reciprocal = alpha.new(alpha_reciprocal.size())
         prev_beta = alpha.new(alpha_reciprocal.size())
 
+    update_tridiag = True
+    last_tridiag_iter = 0
     # Start the iteration
     for k in range(n_iter):
         # Get next alpha
@@ -134,7 +136,7 @@
         # If residual are sufficiently small, then exit loop
         # Alternatively, exit if this is our last iteration
         torch.norm(residual, 2, dim=-2, out=residual_norm)
-        if (residual_norm < tolerance).all() and not n_tridiag:
+        if (residual_norm < tolerance).all() and not (n_tridiag and k < n_tridiag_iter):
             break
 
         # Update precond_residual
@@ -153,7 +155,7 @@
         curr_conjugate_vec.mul_(beta).add_(precond_residual)
 
         # Update tridiagonal matrices, if applicable
-        if n_tridiag and k < n_tridiag_iter:
+        if n_tridiag and k < n_tridiag_iter and update_tridiag:
             alpha_tridiag = alpha.squeeze_(-2).narrow(-1, 0, n_tridiag)
             beta_tridiag = beta.squeeze_(-2).narrow(-1, 0, n_tridiag)
             torch.reciprocal(alpha_tridiag, out=alpha_reciprocal)
@@ -165,6 +167,11 @@
                 torch.mul(prev_beta.sqrt_(), prev_alpha_reciprocal, out=t_mat[k, k - 1])
                 t_mat[k - 1, k].copy_(t_mat[k, k - 1])
 
+                if t_mat[k - 1, k].max() < 1e-6:
+                    update_tridiag = False
+
+            last_tridiag_iter = k
+
             prev_alpha_reciprocal.copy_(alpha_reciprocal)
             prev_beta.copy_(beta_tridiag)
 
@@ -172,6 +179,7 @@
         result = result.squeeze(-1)
 
     if n_tridiag:
+        t_mat = t_mat[:last_tridiag_iter + 1, :last_tridiag_iter + 1]
         if rhs.ndimension() == 3:
             return result, t_mat.permute(3, 2, 0, 1).contiguous()
         else:

{"golden_diff": "diff --git a/gpytorch/utils/linear_cg.py b/gpytorch/utils/linear_cg.py\n--- a/gpytorch/utils/linear_cg.py\n+++ b/gpytorch/utils/linear_cg.py\n@@ -113,6 +113,8 @@\n         prev_alpha_reciprocal = alpha.new(alpha_reciprocal.size())\n         prev_beta = alpha.new(alpha_reciprocal.size())\n \n+    update_tridiag = True\n+    last_tridiag_iter = 0\n     # Start the iteration\n     for k in range(n_iter):\n         # Get next alpha\n@@ -134,7 +136,7 @@\n         # If residual are sufficiently small, then exit loop\n         # Alternatively, exit if this is our last iteration\n         torch.norm(residual, 2, dim=-2, out=residual_norm)\n-        if (residual_norm < tolerance).all() and not n_tridiag:\n+        if (residual_norm < tolerance).all() and not (n_tridiag and k < n_tridiag_iter):\n             break\n \n         # Update precond_residual\n@@ -153,7 +155,7 @@\n         curr_conjugate_vec.mul_(beta).add_(precond_residual)\n \n         # Update tridiagonal matrices, if applicable\n-        if n_tridiag and k < n_tridiag_iter:\n+        if n_tridiag and k < n_tridiag_iter and update_tridiag:\n             alpha_tridiag = alpha.squeeze_(-2).narrow(-1, 0, n_tridiag)\n             beta_tridiag = beta.squeeze_(-2).narrow(-1, 0, n_tridiag)\n             torch.reciprocal(alpha_tridiag, out=alpha_reciprocal)\n@@ -165,6 +167,11 @@\n                 torch.mul(prev_beta.sqrt_(), prev_alpha_reciprocal, out=t_mat[k, k - 1])\n                 t_mat[k - 1, k].copy_(t_mat[k, k - 1])\n \n+                if t_mat[k - 1, k].max() < 1e-6:\n+                    update_tridiag = False\n+\n+            last_tridiag_iter = k\n+\n             prev_alpha_reciprocal.copy_(alpha_reciprocal)\n             prev_beta.copy_(beta_tridiag)\n \n@@ -172,6 +179,7 @@\n         result = result.squeeze(-1)\n \n     if n_tridiag:\n+        t_mat = t_mat[:last_tridiag_iter + 1, :last_tridiag_iter + 1]\n         if rhs.ndimension() == 3:\n             return result, t_mat.permute(3, 2, 0, 1).contiguous()\n         else:\n", "issue": "Exact GP regression Lapack Error syev\nTraining a simple GP regression model produces following runtime error randomly at some seeds: \r\n`RuntimeError: Lapack Error syev : 3 off-diagonal elements didn't converge to zero at /pytorch/aten/src/TH/generic/THTensorLapack.c:388`.\r\nI have been able to reproduce the error with the following code:\r\n\r\n```\r\nimport math\r\nimport torch\r\nimport gpytorch\r\nfrom matplotlib import pyplot as plt\r\n\r\nfrom torch import optim\r\nfrom gpytorch.kernels import RBFKernel\r\nfrom gpytorch.means import ConstantMean\r\nfrom gpytorch.likelihoods import GaussianLikelihood\r\nfrom gpytorch.random_variables import GaussianRandomVariable\r\n\r\ntorch.manual_seed(1)\r\ntrain_x = torch.FloatTensor([[ 6,  6],\r\n       [ 7,  3],\r\n       [ 7,  8],\r\n       [ 8,  6],\r\n       [10,  0]])\r\n\r\ntrain_y = torch.FloatTensor([2.38922077, 0.35414244, 2.24837906, 1.40895557, 0.68164634])\r\n\r\nclass ExactGPModel(gpytorch.models.ExactGP):\r\n    def __init__(self, train_x, train_y, likelihood):\r\n        super(ExactGPModel, self).__init__(train_x, train_y, likelihood)\r\n        self.mean_module = ConstantMean(constant_bounds=(-10, 10))\r\n        self.covar_module = RBFKernel(log_lengthscale_bounds=(-5, 5))\r\n    \r\n    def forward(self, x):\r\n        mean_x = self.mean_module(x)\r\n        covar_x = self.covar_module(x)\r\n        return GaussianRandomVariable(mean_x, covar_x)\r\n\r\nlikelihood = GaussianLikelihood(log_noise_bounds=(-5, 5))\r\nmodel = ExactGPModel(train_x.data, train_y.data, likelihood)\r\n\r\nmodel.train()\r\nlikelihood.train()\r\n\r\noptimizer = torch.optim.Adam([\r\n    {'params': model.parameters()},  \r\n], lr=0.1)\r\n\r\nmll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\r\n\r\ntraining_iter = 50\r\nfor i in range(training_iter):\r\n    optimizer.zero_grad()\r\n    output = model(train_x)\r\n    loss = -mll(output, train_y)\r\n    loss.backward()\r\n    print('Iter %d/%d - Loss: %.3f   log_lengthscale: %.3f   log_noise: %.3f' % (\r\n        i + 1, training_iter, loss.data[0],\r\n        model.covar_module.log_lengthscale.data[0, 0],\r\n        model.likelihood.log_noise.data[0]\r\n    ))\r\n    optimizer.step()\r\n\r\n\r\n```\r\nI am using the latest version of Gpytorch and Pytorch (0.4.0). The above code is identical to the simple gp regression example provided with this repo with changes in train_x and train_y only. Also, the code works fine for seed = 0. \r\n\n", "before_files": [{"content": "from __future__ import absolute_import\nfrom __future__ import division\nfrom __future__ import print_function\nfrom __future__ import unicode_literals\n\nimport torch\nfrom .. import settings\n\n\ndef _default_preconditioner(x):\n    return x.clone()\n\n\ndef linear_cg(\n    matmul_closure,\n    rhs,\n    n_tridiag=0,\n    tolerance=1e-6,\n    eps=1e-20,\n    max_iter=None,\n    max_tridiag_iter=None,\n    initial_guess=None,\n    preconditioner=None,\n):\n    \"\"\"\n    Implements the linear conjugate gradients method for (approximately) solving systems of the form\n\n        lhs result = rhs\n\n    for positive definite and symmetric matrices.\n\n    Args:\n      - matmul_closure - a function which performs a left matrix multiplication with lhs_mat\n      - rhs - the right-hand side of the equation\n      - n_tridiag - returns a tridiagonalization of the first n_tridiag columns of rhs\n      - tolerance - stop the solve when the max residual is less than this\n      - eps - noise to add to prevent division by zero\n      - max_iter - the maximum number of CG iterations\n      - max_tridiag_iter - the maximum size of the tridiagonalization matrix\n      - initial_guess - an initial guess at the solution `result`\n      - precondition_closure - a functions which left-preconditions a supplied vector\n\n    Returns:\n      result - a solution to the system (if n_tridiag is 0)\n      result, tridiags - a solution to the system, and corresponding tridiagonal matrices (if n_tridiag > 0)\n    \"\"\"\n    # Unsqueeze, if necesasry\n    is_vector = rhs.ndimension() == 1\n    if is_vector:\n        rhs = rhs.unsqueeze(-1)\n\n    # Some default arguments\n    if max_iter is None:\n        max_iter = settings.max_cg_iterations.value()\n    if max_tridiag_iter is None:\n        max_tridiag_iter = settings.max_lanczos_quadrature_iterations.value()\n    if initial_guess is None:\n        initial_guess = rhs.new(rhs.size()).zero_()\n    if preconditioner is None:\n        preconditioner = _default_preconditioner\n\n    # If we are running m CG iterations, we obviously can't get more than m Lanczos coefficients\n    if max_tridiag_iter > max_iter:\n        raise RuntimeError(\"Getting a tridiagonalization larger than the number of CG iterations run is not possible!\")\n\n    # Check matmul_closure object\n    if torch.is_tensor(matmul_closure):\n        matmul_closure = matmul_closure.matmul\n    elif not callable(matmul_closure):\n        raise RuntimeError(\"matmul_closure must be a tensor, or a callable object!\")\n\n    # Get some constants\n    n_rows = rhs.size(-2)\n    n_iter = min(max_iter, n_rows)\n    n_tridiag_iter = min(max_tridiag_iter, n_rows)\n\n    # result <- x_{0}\n    result = initial_guess\n\n    # residual: residual_{0} = b_vec - lhs x_{0}\n    residual = rhs - matmul_closure(result)\n\n    # Check for NaNs\n    if not torch.equal(residual, residual):\n        raise RuntimeError(\"NaNs encounterd when trying to perform matrix-vector multiplication\")\n\n    # Sometime we're lucky and the preconditioner solves the system right away\n    residual_norm = residual.norm(2, dim=-2)\n    if (residual_norm < tolerance).all() and not n_tridiag:\n        n_iter = 0  # Skip the iteration!\n\n    # Otherwise, let's define precond_residual and curr_conjugate_vec\n    else:\n        # precon_residual{0} = M^-1 residual_{0}\n        precond_residual = preconditioner(residual)\n        curr_conjugate_vec = precond_residual\n        residual_inner_prod = precond_residual.mul(residual).sum(-2, keepdim=True)\n\n        # Define storage matrices\n        mul_storage = residual.new(residual.size())\n        alpha = residual.new(rhs.size(0), 1, rhs.size(-1)) if rhs.ndimension() == 3 else residual.new(1, rhs.size(-1))\n        beta = alpha.new(alpha.size())\n\n    # Define tridiagonal matrices, if applicable\n    if n_tridiag:\n        if rhs.ndimension() == 3:\n            t_mat = residual.new(n_tridiag_iter, n_tridiag_iter, rhs.size(0), n_tridiag).zero_()\n            alpha_reciprocal = alpha.new(rhs.size(0), n_tridiag)\n        else:\n            t_mat = residual.new(n_tridiag_iter, n_tridiag_iter, n_tridiag).zero_()\n            alpha_reciprocal = alpha.new(n_tridiag)\n\n        prev_alpha_reciprocal = alpha.new(alpha_reciprocal.size())\n        prev_beta = alpha.new(alpha_reciprocal.size())\n\n    # Start the iteration\n    for k in range(n_iter):\n        # Get next alpha\n        # alpha_{k} = (residual_{k-1}^T precon_residual{k-1}) / (p_vec_{k-1}^T mat p_vec_{k-1})\n        mvms = matmul_closure(curr_conjugate_vec)\n        torch.mul(curr_conjugate_vec, mvms, out=mul_storage)\n        torch.sum(mul_storage, -2, keepdim=True, out=alpha)\n        alpha.add_(eps)\n        torch.div(residual_inner_prod, alpha, out=alpha)\n\n        # Update result\n        # result_{k} = result_{k-1} + alpha_{k} p_vec_{k-1}\n        torch.addcmul(result, alpha, curr_conjugate_vec, out=result)\n\n        # Update residual\n        # residual_{k} = residual_{k-1} - alpha_{k} mat p_vec_{k-1}\n        torch.addcmul(residual, -1, alpha, mvms, out=residual)\n\n        # If residual are sufficiently small, then exit loop\n        # Alternatively, exit if this is our last iteration\n        torch.norm(residual, 2, dim=-2, out=residual_norm)\n        if (residual_norm < tolerance).all() and not n_tridiag:\n            break\n\n        # Update precond_residual\n        # precon_residual{k} = M^-1 residual_{k}\n        precond_residual = preconditioner(residual)\n\n        # beta_{k} = (precon_residual{k}^T r_vec_{k}) / (precon_residual{k-1}^T r_vec_{k-1})\n        residual_inner_prod.add_(eps)\n        torch.reciprocal(residual_inner_prod, out=beta)\n        torch.mul(residual, precond_residual, out=mul_storage)\n        torch.sum(mul_storage, -2, keepdim=True, out=residual_inner_prod)\n        beta.mul_(residual_inner_prod)\n\n        # Update curr_conjugate_vec\n        # curr_conjugate_vec_{k} = precon_residual{k} + beta_{k} curr_conjugate_vec_{k-1}\n        curr_conjugate_vec.mul_(beta).add_(precond_residual)\n\n        # Update tridiagonal matrices, if applicable\n        if n_tridiag and k < n_tridiag_iter:\n            alpha_tridiag = alpha.squeeze_(-2).narrow(-1, 0, n_tridiag)\n            beta_tridiag = beta.squeeze_(-2).narrow(-1, 0, n_tridiag)\n            torch.reciprocal(alpha_tridiag, out=alpha_reciprocal)\n\n            if k == 0:\n                t_mat[k, k].copy_(alpha_reciprocal)\n            else:\n                torch.addcmul(alpha_reciprocal, prev_beta, prev_alpha_reciprocal, out=t_mat[k, k])\n                torch.mul(prev_beta.sqrt_(), prev_alpha_reciprocal, out=t_mat[k, k - 1])\n                t_mat[k - 1, k].copy_(t_mat[k, k - 1])\n\n            prev_alpha_reciprocal.copy_(alpha_reciprocal)\n            prev_beta.copy_(beta_tridiag)\n\n    if is_vector:\n        result = result.squeeze(-1)\n\n    if n_tridiag:\n        if rhs.ndimension() == 3:\n            return result, t_mat.permute(3, 2, 0, 1).contiguous()\n        else:\n            return result, t_mat.permute(2, 0, 1).contiguous()\n    else:\n        return result\n", "path": "gpytorch/utils/linear_cg.py"}], "after_files": [{"content": "from __future__ import absolute_import\nfrom __future__ import division\nfrom __future__ import print_function\nfrom __future__ import unicode_literals\n\nimport torch\nfrom .. import settings\n\n\ndef _default_preconditioner(x):\n    return x.clone()\n\n\ndef linear_cg(\n    matmul_closure,\n    rhs,\n    n_tridiag=0,\n    tolerance=1e-6,\n    eps=1e-20,\n    max_iter=None,\n    max_tridiag_iter=None,\n    initial_guess=None,\n    preconditioner=None,\n):\n    \"\"\"\n    Implements the linear conjugate gradients method for (approximately) solving systems of the form\n\n        lhs result = rhs\n\n    for positive definite and symmetric matrices.\n\n    Args:\n      - matmul_closure - a function which performs a left matrix multiplication with lhs_mat\n      - rhs - the right-hand side of the equation\n      - n_tridiag - returns a tridiagonalization of the first n_tridiag columns of rhs\n      - tolerance - stop the solve when the max residual is less than this\n      - eps - noise to add to prevent division by zero\n      - max_iter - the maximum number of CG iterations\n      - max_tridiag_iter - the maximum size of the tridiagonalization matrix\n      - initial_guess - an initial guess at the solution `result`\n      - precondition_closure - a functions which left-preconditions a supplied vector\n\n    Returns:\n      result - a solution to the system (if n_tridiag is 0)\n      result, tridiags - a solution to the system, and corresponding tridiagonal matrices (if n_tridiag > 0)\n    \"\"\"\n    # Unsqueeze, if necesasry\n    is_vector = rhs.ndimension() == 1\n    if is_vector:\n        rhs = rhs.unsqueeze(-1)\n\n    # Some default arguments\n    if max_iter is None:\n        max_iter = settings.max_cg_iterations.value()\n    if max_tridiag_iter is None:\n        max_tridiag_iter = settings.max_lanczos_quadrature_iterations.value()\n    if initial_guess is None:\n        initial_guess = rhs.new(rhs.size()).zero_()\n    if preconditioner is None:\n        preconditioner = _default_preconditioner\n\n    # If we are running m CG iterations, we obviously can't get more than m Lanczos coefficients\n    if max_tridiag_iter > max_iter:\n        raise RuntimeError(\"Getting a tridiagonalization larger than the number of CG iterations run is not possible!\")\n\n    # Check matmul_closure object\n    if torch.is_tensor(matmul_closure):\n        matmul_closure = matmul_closure.matmul\n    elif not callable(matmul_closure):\n        raise RuntimeError(\"matmul_closure must be a tensor, or a callable object!\")\n\n    # Get some constants\n    n_rows = rhs.size(-2)\n    n_iter = min(max_iter, n_rows)\n    n_tridiag_iter = min(max_tridiag_iter, n_rows)\n\n    # result <- x_{0}\n    result = initial_guess\n\n    # residual: residual_{0} = b_vec - lhs x_{0}\n    residual = rhs - matmul_closure(result)\n\n    # Check for NaNs\n    if not torch.equal(residual, residual):\n        raise RuntimeError(\"NaNs encounterd when trying to perform matrix-vector multiplication\")\n\n    # Sometime we're lucky and the preconditioner solves the system right away\n    residual_norm = residual.norm(2, dim=-2)\n    if (residual_norm < tolerance).all() and not n_tridiag:\n        n_iter = 0  # Skip the iteration!\n\n    # Otherwise, let's define precond_residual and curr_conjugate_vec\n    else:\n        # precon_residual{0} = M^-1 residual_{0}\n        precond_residual = preconditioner(residual)\n        curr_conjugate_vec = precond_residual\n        residual_inner_prod = precond_residual.mul(residual).sum(-2, keepdim=True)\n\n        # Define storage matrices\n        mul_storage = residual.new(residual.size())\n        alpha = residual.new(rhs.size(0), 1, rhs.size(-1)) if rhs.ndimension() == 3 else residual.new(1, rhs.size(-1))\n        beta = alpha.new(alpha.size())\n\n    # Define tridiagonal matrices, if applicable\n    if n_tridiag:\n        if rhs.ndimension() == 3:\n            t_mat = residual.new(n_tridiag_iter, n_tridiag_iter, rhs.size(0), n_tridiag).zero_()\n            alpha_reciprocal = alpha.new(rhs.size(0), n_tridiag)\n        else:\n            t_mat = residual.new(n_tridiag_iter, n_tridiag_iter, n_tridiag).zero_()\n            alpha_reciprocal = alpha.new(n_tridiag)\n\n        prev_alpha_reciprocal = alpha.new(alpha_reciprocal.size())\n        prev_beta = alpha.new(alpha_reciprocal.size())\n\n    update_tridiag = True\n    last_tridiag_iter = 0\n    # Start the iteration\n    for k in range(n_iter):\n        # Get next alpha\n        # alpha_{k} = (residual_{k-1}^T precon_residual{k-1}) / (p_vec_{k-1}^T mat p_vec_{k-1})\n        mvms = matmul_closure(curr_conjugate_vec)\n        torch.mul(curr_conjugate_vec, mvms, out=mul_storage)\n        torch.sum(mul_storage, -2, keepdim=True, out=alpha)\n        alpha.add_(eps)\n        torch.div(residual_inner_prod, alpha, out=alpha)\n\n        # Update result\n        # result_{k} = result_{k-1} + alpha_{k} p_vec_{k-1}\n        torch.addcmul(result, alpha, curr_conjugate_vec, out=result)\n\n        # Update residual\n        # residual_{k} = residual_{k-1} - alpha_{k} mat p_vec_{k-1}\n        torch.addcmul(residual, -1, alpha, mvms, out=residual)\n\n        # If residual are sufficiently small, then exit loop\n        # Alternatively, exit if this is our last iteration\n        torch.norm(residual, 2, dim=-2, out=residual_norm)\n        if (residual_norm < tolerance).all() and not (n_tridiag and k < n_tridiag_iter):\n            break\n\n        # Update precond_residual\n        # precon_residual{k} = M^-1 residual_{k}\n        precond_residual = preconditioner(residual)\n\n        # beta_{k} = (precon_residual{k}^T r_vec_{k}) / (precon_residual{k-1}^T r_vec_{k-1})\n        residual_inner_prod.add_(eps)\n        torch.reciprocal(residual_inner_prod, out=beta)\n        torch.mul(residual, precond_residual, out=mul_storage)\n        torch.sum(mul_storage, -2, keepdim=True, out=residual_inner_prod)\n        beta.mul_(residual_inner_prod)\n\n        # Update curr_conjugate_vec\n        # curr_conjugate_vec_{k} = precon_residual{k} + beta_{k} curr_conjugate_vec_{k-1}\n        curr_conjugate_vec.mul_(beta).add_(precond_residual)\n\n        # Update tridiagonal matrices, if applicable\n        if n_tridiag and k < n_tridiag_iter and update_tridiag:\n            alpha_tridiag = alpha.squeeze_(-2).narrow(-1, 0, n_tridiag)\n            beta_tridiag = beta.squeeze_(-2).narrow(-1, 0, n_tridiag)\n            torch.reciprocal(alpha_tridiag, out=alpha_reciprocal)\n\n            if k == 0:\n                t_mat[k, k].copy_(alpha_reciprocal)\n            else:\n                torch.addcmul(alpha_reciprocal, prev_beta, prev_alpha_reciprocal, out=t_mat[k, k])\n                torch.mul(prev_beta.sqrt_(), prev_alpha_reciprocal, out=t_mat[k, k - 1])\n                t_mat[k - 1, k].copy_(t_mat[k, k - 1])\n\n                if t_mat[k - 1, k].max() < 1e-6:\n                    update_tridiag = False\n\n            last_tridiag_iter = k\n\n            prev_alpha_reciprocal.copy_(alpha_reciprocal)\n            prev_beta.copy_(beta_tridiag)\n\n    if is_vector:\n        result = result.squeeze(-1)\n\n    if n_tridiag:\n        t_mat = t_mat[:last_tridiag_iter + 1, :last_tridiag_iter + 1]\n        if rhs.ndimension() == 3:\n            return result, t_mat.permute(3, 2, 0, 1).contiguous()\n        else:\n            return result, t_mat.permute(2, 0, 1).contiguous()\n    else:\n        return result\n", "path": "gpytorch/utils/linear_cg.py"}]}