def solve_G(self, b):
"""Solve ``Gx=b`` where ``G`` is the Gram matrix of the GP.
Uses the internally-stored LU decomposition of ``G`` computed in
``gp.update()``."""
assert self.ready
return linalg.lu_solve((self.LU, self.LU_piv), b, check_finite=True)
python类lu_solve()的实例源码
gaussian_process.py 文件源码
项目:probabilistic_line_search
作者: ProbabilisticNumerics
项目源码
文件源码
阅读 38
收藏 0
点赞 0
评论 0
def lu_solve_ATAI(A, rho, b, lu, piv):
r"""
Solve the linear system :math:`(A^T A + \rho I)\mathbf{x} = \mathbf{b}`
or :math:`(A^T A + \rho I)X = B` using :func:`scipy.linalg.lu_solve`.
Parameters
----------
A : array_like
Matrix :math:`A`
rho : float
Scalar :math:`\rho`
b : array_like
Vector :math:`\mathbf{b}` or matrix :math:`B`
lu : array_like
Matrix containing U in its upper triangle, and L in its lower triangle,
as returned by :func:`scipy.linalg.lu_factor`
piv : array_like
Pivot indices representing the permutation matrix P, as returned by
:func:`scipy.linalg.lu_factor`
Returns
-------
x : ndarray
Solution to the linear system.
"""
N, M = A.shape
if N >= M:
x = linalg.lu_solve((lu, piv), b)
else:
x = (b - A.T.dot(linalg.lu_solve((lu, piv), A.dot(b), 1))) / rho
return x
def lu_solve_AATI(A, rho, b, lu, piv):
r"""
Solve the linear system :math:`(A A^T + \rho I)\mathbf{x} = \mathbf{b}`
or :math:`(A A^T + \rho I)X = B` using :func:`scipy.linalg.lu_solve`.
Parameters
----------
A : array_like
Matrix :math:`A`
rho : float
Scalar :math:`\rho`
b : array_like
Vector :math:`\mathbf{b}` or matrix :math:`B`
lu : array_like
Matrix containing U in its upper triangle, and L in its lower triangle,
as returned by :func:`scipy.linalg.lu_factor`
piv : array_like
Pivot indices representing the permutation matrix P, as returned by
:func:`scipy.linalg.lu_factor`
Returns
-------
x : ndarray
Solution to the linear system.
"""
N, M = A.shape
if N >= M:
x = (b - linalg.lu_solve((lu, piv), b.dot(A).T).T.dot(A.T)) / rho
else:
x = linalg.lu_solve((lu, piv), b.T).T
return x
def solve(self, rhs):
"""Solve with the system matrix
Args:
rhs: Right hand side in system
Returns:
solution to M*x = rhs
"""
self.update()
return la.lu_solve(self.lupiv, rhs)
def mysolve(LU, b):
if isinstance(LU, SuperLU):
return LU.solve(b)
else:
return lu_solve(LU, b)
###############################################################################
# Solve via full system
###############################################################################
def solve_nonlinear(self, inputs, outputs):
"""
Use numpy to solve Ax=b for x.
Parameters
----------
inputs : Vector
unscaled, dimensional input variables read via inputs[key]
outputs : Vector
unscaled, dimensional output variables read via outputs[key]
"""
# lu factorization for use with solve_linear
self._lup = linalg.lu_factor(inputs['A'])
outputs['x'] = linalg.lu_solve(self._lup, inputs['b'])
def solve_linear(self, d_outputs, d_residuals, mode):
r"""
Back-substitution to solve the derivatives of the linear system.
If mode is:
'fwd': d_residuals \|-> d_outputs
'rev': d_outputs \|-> d_residuals
Parameters
----------
d_outputs : Vector
unscaled, dimensional quantities read via d_outputs[key]
d_residuals : Vector
unscaled, dimensional quantities read via d_residuals[key]
mode : str
either 'fwd' or 'rev'
"""
if mode == 'fwd':
sol_vec, rhs_vec = d_outputs, d_residuals
t = 0
else:
sol_vec, rhs_vec = d_residuals, d_outputs
t = 1
# print("foobar", rhs_vec['x'])
sol_vec['x'] = linalg.lu_solve(self._lup, rhs_vec['x'], trans=t)
def solve_nonlinear(self, inputs, outputs):
force_vector = np.concatenate([self.metadata['force_vector'], np.zeros(2)])
self.lu = lu_factor(inputs['K'])
outputs['d'] = lu_solve(self.lu, force_vector)
def solve_linear(self, d_outputs, d_residuals, mode):
if mode == 'fwd':
d_outputs['d'] = lu_solve(self.lu, d_residuals['d'], trans=0)
else:
d_residuals['d'] = lu_solve(self.lu, d_outputs['d'], trans=1)
def _matvec(self, x):
return lu_solve(self.M_lu, x)