import logging
import numpy as np
__docformat__ = 'restructuredtext'
# Default (null) logger.
null_log = logging.getLogger('krylov')
null_log.setLevel(logging.INFO)
null_log.addHandler(logging.NullHandler())
class KrylovMethod(object):
def __init__(self, op, **kwargs):
self.prefix = 'Generic: '
self.name = 'Generic Krylov Method (must be subclassed)'
# Mandatory arguments
self.op = op
# Optional keyword arguments
self.abstol = kwargs.get('abstol', 1.0e-8)
self.reltol = kwargs.get('reltol', 1.0e-6)
self.precon = kwargs.get('precon', None)
self.logger = kwargs.get('logger', null_log)
self.residNorm = None
self.residNorm0 = None
self.residHistory = []
self.nMatvec = 0
self.nIter = 0
self.converged = False
self.bestSolution = None
self.x = self.bestSolution
def _write(self, msg):
# If levels other than info are needed they should be used explicitly.
self.logger.info(msg)
def solve(self, rhs, **kwargs):
"""
This is the :meth:`solve` method of the abstract `KrylovMethod` class.
The class must be specialized and this method overridden.
"""
# raise NotImplementedError('This method must be subclassed')
"""
A Python implementation of SYMMLQ.
This is a line-by-line translation from Matlab code
available at http://[Log in to view URL]
.. moduleauthor: D. Orban <dominique.orban@gerad.ca>
"""
__docformat__ = 'restructuredtext'
def machine_epsilon():
"""Return the machine epsilon in double precision."""
return np.finfo(np.double).eps
class Symmlq(KrylovMethod) :
def __init__(self, op, **kwargs):
KrylovMethod.__init__(self, op, **kwargs)
self.name = 'Symmetric Indefinite Lanczos with Orthogonal Factorization'
self.acronym = 'SYMMLQ'
self.prefix = self.acronym + ': '
self.iterates = []
def solve(self, rhs, **kwargs):
# Set parameters
n = rhs.shape[0]
nMatvec = 0
matvec_max = kwargs.get('matvec_max', 2*n+2)
rtol = kwargs.get('rtol', 1.0e-9)
check = kwargs.get('check', False)
shift = kwargs.get('shift', None)
if shift == 0.0: shift = None
eps = machine_epsilon()
store_iterates = kwargs.get('store_iterates', False)
first = 'Enter SYMMLQ. '
last = 'Exit SYMMLQ. '
space = ' '
msg={
-1:' beta2 = 0. If M = I, b and x are eigenvectors',
0:' beta1 = 0. The exact solution is x = 0',
1:' Requested accuracy achieved, as determined by rtol',
2:' Reasonable accuracy achieved, given eps',
3:' x has converged to an eigenvector',
4:' acond has exceeded 0.1/eps',
5:' The iteration limit was reached',
6:' aprod does not define a symmetric matrix',
7:' msolve does not define a symmetric matrix',
8:' msolve does not define a pos-def preconditioner'}
self.logger.info(first + 'Solution of symmetric Ax = b')
fmt = 'n = %3g precon = %5s '
self.logger.info(fmt % (n, repr(self.precon is None)))
if shift is not None: self.logger.info('shift = %23.14e' % shift)
fmt = 'maxit = %3g eps = %11.2e rtol = %11.2e'
self.logger.info(fmt % (int((matvec_max-2.0)/2),eps,rtol))
istop = 0 ; ynorm = 0 ; w = np.zeros(n) ; acond = 0
itn = 0 ; xnorm = 0 ; x = np.zeros(n) ; done=False
anorm = 0 ; rnorm = 0 ; v = np.zeros(n)
if store_iterates:
self.iterates.append(x.copy())
# Set up y for the first Lanczos vector v1.
# y is really beta1 * P * v1 where P = C^(-1).
# y and beta1 will be zero if b = 0.
r1 = rhs.copy()
if self.precon is not None:
y = self.precon * r1
else:
y = rhs.copy()
b1 = y[0] ; beta1 = np.dot(r1, y)
# Ensure preconditioner is symmetric.
if check and self.precon is not None:
r2 = self.precon * y
s = np.dot(y,y)
t = np.dot(r1,r2)
z = np.abs(s-t)
epsa = (s+eps) * eps**(1.0/3)
if z > epsa:
istop = 7
done = True
# Test for an indefinite preconditioner.
# If rhs = 0 exactly, stop with x = 0.
if beta1 < 0:
istop = 8
done = True
if beta1 == 0:
done = True
if beta1 > 0:
beta1 = np.sqrt(beta1)
s = 1.0 / beta1
v = s * y
y = self.op(v) ; nMatvec += 1
if check:
r2 = self.op * y # Do not count this matrix-vector product
s = np.dot(y,y)
t = np.dot(v,r2)
z = abs(s-t)
epsa = (s+eps) * eps**(1.0/3)
if z > epsa:
istop = 6
done = True
# Set up y for the second Lanczos vector.
# Again, y is beta * P * v2 where P = C^(-1).
# y and beta will be zero or very small if Abar = I or constant * I.
if shift is not None: y -= shift * v
alfa = np.dot(v, y)
y -= (alfa / beta1) * r1
# Make sure r2 will be orthogonal to the first v.
z = np.dot(v, y)
s = np.dot(v, v)
y -= (z / s) * v
r2 = y.copy()
if self.precon is not None: y = self.precon * r2
oldb = beta1
beta = np.dot(r2, y)
if beta < 0:
istop = 8
done = True
# Cause termination (later) if beta is essentially zero.
beta = np.sqrt(beta)
if beta <= eps:
istop = -1
# See if the local reorthogonalization achieved anything.
denom = np.sqrt(s) * np.linalg.norm(r2) + eps
s = z / denom
t = np.dot(v, r2)
t = t / denom
self.logger.info('beta1 = %10.2e alpha1 = %9.2e'% (beta1,alfa))
self.logger.info('(v1, v2) before and after %14.2e' % s)
self.logger.info('local reorthogonalization %14.2e' % t)
# Initialize other quantities.
cgnorm = beta1 ; rhs2 = 0 ; tnorm = alfa**2 + beta**2
gbar = alfa ; bstep = 0 ; ynorm2 = 0
dbar = beta ; snprod = 1 ; gmax = np.abs(alfa) + eps
rhs1 = beta1 ; x1cg = 0 ; gmin = gmax
qrnorm = beta1
# end if beta1 > 0
head1 = ' Itn x(1)(cg) normr(cg) r(minres)'
head2 = ' bstep anorm acond'
self.logger.info(head1 + head2)
str1 = '%6g %12.5e %10.3e' % (itn, x1cg, cgnorm)
str2 = ' %10.3e %8.1e' % (qrnorm, bstep/beta1)
self.logger.info(str1 + str2)
# ------------------------------------------------------------------
# Main iteration loop.
# ------------------------------------------------------------------
# Estimate various norms and test for convergence.
if not done:
while nMatvec < matvec_max:
itn = itn + 1
anorm = np.sqrt(tnorm)
ynorm = np.sqrt(ynorm2)
epsa = anorm * eps
epsx = anorm * ynorm * eps
epsr = anorm * ynorm * rtol
diag = gbar
if diag == 0: diag = epsa
lqnorm = np.sqrt(rhs1**2 + rhs2**2)
qrnorm = snprod * beta1
cgnorm = qrnorm * beta / np.abs(diag)
# Estimate Cond(A).
# In this version we look at the diagonals of L in the
# factorization of the tridiagonal matrix, T = L*Q.
# Sometimes, T(k) can be misleadingly ill-conditioned when
# T(k+1) is not, so we must be careful not to overestimate acond
if lqnorm < cgnorm:
acond = gmax / gmin
else:
denom = min(gmin, np.abs(diag))
acond = gmax / denom
zbar = rhs1 / diag
z = (snprod * zbar + bstep) / beta1
x1lq = x[0] + b1 * bstep / beta1
x1cg = x[0] + w[0] * zbar + b1 * z
# See if any of the stopping criteria are satisfied.
# In rare cases, istop is already -1 from above
# (Abar = const * I).
if istop == 0:
if nMatvec >= matvec_max : istop = 5
if acond >= 0.1/eps : istop = 4
if epsx >= beta1 : istop = 3
if cgnorm <= epsx : istop = 2
if cgnorm <= epsr : istop = 1
prnt = False
if n <= 40 : prnt = True
if nMatvec <= 20 : prnt = True
if nMatvec >= matvec_max - 10 : prnt = True
if itn%10 == 0 : prnt = True
if cgnorm <= 10.0*epsx : prnt = True
if cgnorm <= 10.0*epsr : prnt = True
if acond >= 0.01/eps : prnt = True
if istop != 0 : prnt = True
if prnt:
str1 = '%6g %12.5e %10.3e' % (itn, x1cg, cgnorm)
str2 = ' %10.3e %8.1e' % (qrnorm, bstep/beta1)
str3 = ' %8.1e %8.1e' % (anorm, acond)
self.logger.info(str1 + str2 + str3)
if istop !=0:
break
# Obtain the current Lanczos vector v = (1 / beta)*y
# and set up y for the next iteration.
s = 1/beta
v = s * y
y = self.op (v) ; nMatvec += 1
if shift is not None: y -= shift * v
y -= (beta / oldb) * r1
alfa = np.dot(v, y)
y -= (alfa / beta) * r2
r1 = r2.copy()
r2 = y.copy()
if self.precon is not None: y = self.precon * r2
oldb = beta
beta = np.dot(r2, y)
if beta < 0:
istop = 6
break
beta = np.sqrt(beta)
tnorm = tnorm + alfa**2 + oldb**2 + beta**2
# Compute the next plane rotation for Q.
gamma = np.sqrt(gbar**2 + oldb**2)
cs = gbar / gamma
sn = oldb / gamma
delta = cs * dbar + sn * alfa
gbar = sn * dbar - cs * alfa
epsln = sn * beta
dbar = - cs * beta
# Update X.
z = rhs1 / gamma
s = z*cs
t = z*sn
x += s*w + t*v
w *= sn ; w -= cs*v
if store_iterates:
self.iterates.append(x.copy())
# Accumulate the step along the direction b, and go round again.
bstep = snprod * cs * z + bstep
snprod = snprod * sn
gmax = max(gmax, gamma)
gmin = min(gmin, gamma)
ynorm2 = z**2 + ynorm2
rhs1 = rhs2 - delta * z
rhs2 = - epsln * z
# end while
# ------------------------------------------------------------------
# End of main iteration loop.
# ------------------------------------------------------------------
# Move to the CG point if it seems better.
# In this version of SYMMLQ, the convergence tests involve
# only cgnorm, so we're unlikely to stop at an LQ point,
# EXCEPT if the iteration limit interferes.
if cgnorm < lqnorm:
zbar = rhs1 / diag
bstep = snprod * zbar + bstep
ynorm = np.sqrt(ynorm2 + zbar**2)
x += zbar * w
# Add the step along b.
bstep = bstep / beta1
if self.precon is not None:
y = self.precon * rhs
else:
y = rhs.copy()
x += bstep * y
return x
def sparse_prod(spA, x):
y = []
for i in range(len(x)):
xi = sum(x[j] * v for j, v in spA[i].items())
y.append(xi)
return y
A_dct = {
0: {0: 2, 1: 3, 2: 5},
1: {0: 3, 1: 2, 2: 1},
2: {0: 5, 1: 1, 2: 4},
}
A = np.array([[2,3,5],[3,2,1],[5,1,4]])
b = np.array([1,2,3])
kr = Symmlq(lambda x: sparse_prod(A_dct, x))
ans = kr.solve(b)
print(ans)
print(np.linalg.solve(A, b))
To embed this program on your website, copy the following code and paste it into your website's HTML: