python类lu_factor()的实例源码

def update(self):
    """Set up the Gram matrix and compute its LU decomposition to make the GP
    ready for inference (calls to ``.gp.mu(t)``, ``gp.V(t)``, etc...).

    Call this method after you have manipulated the GP by
       - ``gp.reset()`` ing,
       - adding observations with ``gp.add(t, f, df)``, or
       - adjusting the sigmas via ``gp.update_sigmas()``.
    and want to perform inference next."""

    if self.ready:
      return

    # Set up the kernel matrices.
    self.K = np.matrix(np.zeros([self.N, self.N]))
    self.Kd = np.matrix(np.zeros([self.N, self.N]))
    self.dKd = np.matrix(np.zeros([self.N, self.N]))    
    for i in range(self.N):
      for j in range(self.N):
        self.K[i, j] = self.k(self.ts[i], self.ts[j])
        self.Kd[i, j] = self.kd(self.ts[i], self.ts[j])
        self.dKd[i, j] = self.dkd(self.ts[i], self.ts[j])

    # Put together the Gram matrix
    S_f = np.matrix(np.diag(self.fvars))
    S_df = np.matrix(np.diag(self.dfvars))
    self.G = np.bmat([[self.K + S_f, self.Kd],
                      [self.Kd.T, self.dKd + S_df]])

    # Compute the LU decomposition of G and store it
    self.LU, self.LU_piv = linalg.lu_factor(self.G, check_finite=True)

    # Set ready switch to True
    self.ready = True

    # Pre-compute the regression weights used in mu
    self.w = self.solve_G(np.array(self.fs + self.dfs))

def lu_factor(A, rho):
    r"""
    Compute LU factorisation of either :math:`A^T A + \rho I` or
    :math:`A A^T + \rho I`, depending on which matrix is smaller.

    Parameters
    ----------
    A : array_like
      Array :math:`A`
    rho : float
      Scalar :math:`\rho`

    Returns
    -------
    lu : ndarray
      Matrix containing U in its upper triangle, and L in its lower triangle,
      as returned by :func:`scipy.linalg.lu_factor`
    piv : ndarray
      Pivot indices representing the permutation matrix P, as returned by
      :func:`scipy.linalg.lu_factor`
    """

    N, M = A.shape
    # If N < M it is cheaper to factorise A*A^T + rho*I and then use the
    # matrix inversion lemma to compute the inverse of A^T*A + rho*I
    if N >= M:
        lu, piv = linalg.lu_factor(A.T.dot(A) + rho*np.identity(M,
                                   dtype=A.dtype))
    else:
        lu, piv = linalg.lu_factor(A.dot(A.T) + rho*np.identity(N,
                                   dtype=A.dtype))
    return lu, piv

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 refactor(self):
        """Compute factorization
        """
        eta = min(1e-5, 1e-16 * np.sqrt(la.norm(self.M, 1) * la.norm(self.M, np.inf))) * np.eye(self.M.shape[0])
        self.lupiv = la.lu_factor(self.M + eta)

def myfactor(A):
    if issparse(A):
        return splu(A.tocsc())
    else:
        return lu_factor(A)

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_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 linearize(self, inputs, outputs, partials):
        num_elements = self.metadata['num_elements']
        num_nodes = num_elements + 1
        size = 2 * num_nodes + 2

        self.lu = lu_factor(inputs['K'])

        partials['d', 'K'] = np.outer(np.ones(size), outputs['d']).flatten()
        partials['d', 'd'] = inputs['K']

def __init__(self, M):
        self.M_lu = lu_factor(M)
        self.shape = M.shape
        self.dtype = M.dtype