Source code for mmf.math.linalg

"""Linear Algebra Routines"""
from __future__ import division

__all__ = ['modchol', 'paradiso', 'block_diag']

import numpy as np
from cholesky.gmw81 import colmod as modchol

def block_diag(arrays):
    """Create a new diagonal matrix from the provided arrays.

    Parameters
    ----------
    a, b, c, ... : ndarray
        Input arrays.

    Returns
    -------
    D : ndarray
        Array with a, b, c, ... on the diagonal.

    Examples
    --------
    >>> 
    """
    arrays = [np.asarray(a) for a in arrays]
    shapes = np.array([a.shape for a in arrays])
    out = np.zeros(np.sum(shapes, axis=0))

    r, c = 0, 0
    for i, (rr, cc) in enumerate(shapes):
        out[r:r + rr, c:c + cc] = arrays[i]
        r += rr
        c += cc
    return out

def signm(X, accel_tol=1e-2, maxit=16):
    r"""Uses Newton's method to compute the sign function of a matrix.

    From Higham's Matrix Function Toolbox.

    Note: this is only useful for non-Hermitian matrices.  It is generally
    faster to compute the spectrum.
    """
    tol = math.sqrt(len(X))*_EPS/2.0
    need_accel = True
    reldiff = np.inf
    for _k in xrange(maxit):
        Xold = X
        Xinv = np.linalg.inv(X)
        X = 0.5*(X + Xinv)
        diff_F = np.linalg.norm(X - Xold, 'fro')
        reldiff_old = reldiff
        reldiff = diff_F/np.linalg.norm(X,'fro')
        if need_accel and (reldiff < accel_tol):
            need_accel = False
        if (diff_F <= math.sqrt(
                tol*np.linalg.norm(X,'fro')/np.linalg.norm(Xinv,'fro'))
            or
            (reldiff > reldiff_old/2.0 and not need_accel)):
            return X
    raise ValueError("Failed to converge after %i iterations" % (_k,))