Python linalg.hilbert函数代码示例

    def test_basic(self):
        h3 = array([[1.0, 1 / 2.0, 1 / 3.0], [1 / 2.0, 1 / 3.0, 1 / 4.0], [1 / 3.0, 1 / 4.0, 1 / 5.0]])
        assert_array_almost_equal(hilbert(3), h3)

        assert_array_equal(hilbert(1), [[1.0]])

        h0 = hilbert(0)
        assert_equal(h0.shape, (0, 0))

 def test_inverse(self):
     for n in xrange(1, 10):
         a = hilbert(n)
         b = invhilbert(n)
         # The Hilbert matrix is increasingly badly conditioned,
         # so take that into account in the test
         c = cond(a)
         assert_allclose(a.dot(b), eye(n), atol=1e-15*c, rtol=1e-15*c)

def test_svd_which():
    # check that the which parameter works as expected
    x = hilbert(6)
    for which in ['LM', 'SM']:
        _, s, _ = sorted_svd(x, 2, which=which)
        ss = svds(x, 2, which=which, return_singular_vectors=False)
        ss.sort()
        assert_allclose(s, ss, atol=np.sqrt(1e-15))

def test_svd_maxiter():
    # check that maxiter works as expected
    x = hilbert(6)
    # ARPACK shouldn't converge on such an ill-conditioned matrix with just
    # one iteration
    assert_raises(ArpackNoConvergence, svds, x, 1, maxiter=1)
    # but 100 iterations should be more than enough
    u, s, vt = svds(x, 1, maxiter=100)
    assert_allclose(s, [1.7], atol=0.5)

def run_stats(dimension):
    matrices = {}
    # Generation of the 4 relevant matrices
    matrices['normal'] = np.matrix(np.random.randn(dimension, dimension))
    matrices['hilbert'] = np.matrix(sp.hilbert(dimension))
    matrices['pascal'] = np.matrix(pascal(dimension))
    if dimension != 100:
        matrices['magic'] = np.matrix(magic(dimension))
    else:
        matrices['magic'] = np.matrix(np.identity(100))

    x = np.matrix(np.ones((dimension,1)))

    data = {} # these are just bookkeeping things that keep track of variables for future use
    for name, matr in matrices.items():
        this_matrix_data = {}

        b = matr * x
        try:
            xhat = sp.solve(matr, b)
            xhat1 = np.linalg.inv(matr) * b
        except np.linalg.LinAlgError:
            xhat = x
            xhat1 = x

        this_matrix_data['delta_b'] = matr * xhat - b
        norm_data = {}

        # Loop through the three norms we're asked to do
        for norm_type in [1, 2, np.inf]:
            values = {}
            values['x_relative_error'] = sp.norm(x - xhat, norm_type) / sp.norm(x, norm_type)
            values['x_relative_error_inv'] = sp.norm(x - xhat1, norm_type) / sp.norm(x, norm_type)

            try:
                values['condition_no'] = np.linalg.cond(matr, norm_type)
            except:
                values['condition_no'] = 0

            values['cond_rel_b_err'] = values['condition_no'] * (sp.norm(this_matrix_data['delta_b'], norm_type) / sp.norm(b, norm_type))
            norm_data[norm_type] = values

        this_matrix_data['norm_dep_vals'] = norm_data
        data[name] = this_matrix_data

    return data

def demo(n):
    from numpy import random
    from scipy import linalg
    A = random.randn(n, n)
    x = random.randn(n)
    
    print('\n\t\t Random matrix of dimension', n)
    showerr(A, x)

    A = linalg.hilbert(n)
    print('\n\t\t Hilbert matrix of dimension', n)
    showerr(A, x)

    from . import cond
    A = cond.laplace(n)
    print('\n\t\t Disc Laplacian of dimension', n)
    showerr(A, x)

def main():
    for n in range(12, 20):
        H = hilbert(n)
        x = [1000 for _ in range(n)]
        b = H.dot(x)

        x = np.linalg.solve(H, b)

        Hx = H.dot(x)
        r = b - Hx

        print np.linalg.norm(r, np.inf)

        iter = 0
        while np.linalg.norm(r, np.inf) > 1e-5 and iter < 50:
            z = np.linalg.solve(H, r)
            x = x + z
            Ax = H.dot(x)
            r = b - Ax
            iter += 1

        print iter, x

print(A)


#b)
b=np.array([npr.rand(),npr.rand(),npr.rand()])
print(b)
#c)
x=solve(A,b)
print(x)

#sprawdzenie poprawności rozwiązania
print(npl.norm(np.dot(A,x)-b))
# powinno wyjść coś bliskiego zeru

#d)
H=sclin.hilbert(4)
print(H)

#e)
c=np.array([npr.rand(),npr.rand(),npr.rand(),npr.rand()])
print(c)

#f)
y=solve(H,c)
print(y)

#sprawdzenie poprawności rozwiązania
print(npl.norm(np.dot(H,y)-c))
# powinno wyjść coś bliskiego zeru

#g

def test_svd_return():
    # check that the return_singular_vectors parameter works as expected
    x = hilbert(6)
    _, s, _ = sorted_svd(x, 2)
    ss = svds(x, 2, return_singular_vectors=False)
    assert_allclose(s, ss)

B


# Out[3]:

#     array([[ 1.,  2.,  3.,  4.],
#            [ 2.,  3.,  4.,  5.],
#            [ 3.,  4.,  5.,  6.],
#            [ 4.,  5.,  6.,  7.]])

# In[4]:

A = 1/B

from scipy.linalg import hilbert
A2 = hilbert(n)

print "difference: %g" % la.norm(A-A2)
print "cond: %g" % np.linalg.cond(A)


# Out[4]:

#     difference: 0
#     cond: 15513.7
# 

#     /usr/lib/pymodules/python2.7/numpy/oldnumeric/__init__.py:11: ModuleDeprecationWarning: The oldnumeric module will be dropped in Numpy 1.9
#       warnings.warn(_msg, ModuleDeprecationWarning)
# 

    for k in range(n-1,-1,-1):
        b[k] = (b[k] - np.dot(L[k+1:n,k],b[k+1:n]))/L[k,k]
    return b

choleski(a2)
print ("The answer for problem 11 is:")
print (choleskiSol(a2,b2))
"""





#Problem 15
#Code is not working correctly. I did not have enough time to finish it.
a3 = hilbert(2) 
b3 = []
x = []
i = 0
actNorm = 1*10**-6

while True:
    b3.append(a3[i].sum())
    x.append([1])
    xNorm = np.linalg.norm(x, np.inf)
    xApprox = np.linalg.solve(a3, b3)
    xHat = np.linalg.norm(xApprox, np.inf)
    approxErr = xNorm - xHat
    if abs(approxErr) < actNorm:
        break
    i = i + 1

 def test_badcall(self):
     A = hilbert(5).astype(np.float32)
     assert_raises(ValueError, pymatrixid.interp_decomp, A, 1e-6, rand=False)

    def check_id(self, dtype):
        # Test ID routines on a Hilbert matrix.

        # set parameters
        n = 300
        eps = 1e-12

        # construct Hilbert matrix
        A = hilbert(n).astype(dtype)
        if np.issubdtype(dtype, np.complexfloating):
            A = A * (1 + 1j)
        L = aslinearoperator(A)

        # find rank
        S = np.linalg.svd(A, compute_uv=False)
        try:
            rank = np.nonzero(S < eps)[0][0]
        except:
            rank = n

        # print input summary
        _debug_print("Hilbert matrix dimension:        %8i" % n)
        _debug_print("Working precision:               %8.2e" % eps)
        _debug_print("Rank to working precision:       %8i" % rank)

        # set print format
        fmt = "%8.2e (s) / %5s"

        # test real ID routines
        _debug_print("-----------------------------------------")
        _debug_print("Real ID routines")
        _debug_print("-----------------------------------------")

        # fixed precision
        _debug_print("Calling iddp_id / idzp_id  ...",)
        t0 = time.clock()
        k, idx, proj = pymatrixid.interp_decomp(A, eps, rand=False)
        t = time.clock() - t0
        B = pymatrixid.reconstruct_matrix_from_id(A[:, idx[:k]], idx, proj)
        _debug_print(fmt % (t, np.allclose(A, B, eps)))
        assert_(np.allclose(A, B, eps))

        _debug_print("Calling iddp_aid / idzp_aid ...",)
        t0 = time.clock()
        k, idx, proj = pymatrixid.interp_decomp(A, eps)
        t = time.clock() - t0
        B = pymatrixid.reconstruct_matrix_from_id(A[:, idx[:k]], idx, proj)
        _debug_print(fmt % (t, np.allclose(A, B, eps)))
        assert_(np.allclose(A, B, eps))

        _debug_print("Calling iddp_rid / idzp_rid ...",)
        t0 = time.clock()
        k, idx, proj = pymatrixid.interp_decomp(L, eps)
        t = time.clock() - t0
        B = pymatrixid.reconstruct_matrix_from_id(A[:, idx[:k]], idx, proj)
        _debug_print(fmt % (t, np.allclose(A, B, eps)))
        assert_(np.allclose(A, B, eps))

        # fixed rank
        k = rank

        _debug_print("Calling iddr_id / idzr_id  ...",)
        t0 = time.clock()
        idx, proj = pymatrixid.interp_decomp(A, k, rand=False)
        t = time.clock() - t0
        B = pymatrixid.reconstruct_matrix_from_id(A[:, idx[:k]], idx, proj)
        _debug_print(fmt % (t, np.allclose(A, B, eps)))
        assert_(np.allclose(A, B, eps))

        _debug_print("Calling iddr_aid / idzr_aid ...",)
        t0 = time.clock()
        idx, proj = pymatrixid.interp_decomp(A, k)
        t = time.clock() - t0
        B = pymatrixid.reconstruct_matrix_from_id(A[:, idx[:k]], idx, proj)
        _debug_print(fmt % (t, np.allclose(A, B, eps)))
        assert_(np.allclose(A, B, eps))

        _debug_print("Calling iddr_rid / idzr_rid ...",)
        t0 = time.clock()
        idx, proj = pymatrixid.interp_decomp(L, k)
        t = time.clock() - t0
        B = pymatrixid.reconstruct_matrix_from_id(A[:, idx[:k]], idx, proj)
        _debug_print(fmt % (t, np.allclose(A, B, eps)))
        assert_(np.allclose(A, B, eps))

        # check skeleton and interpolation matrices
        idx, proj = pymatrixid.interp_decomp(A, k, rand=False)
        P = pymatrixid.reconstruct_interp_matrix(idx, proj)
        B = pymatrixid.reconstruct_skel_matrix(A, k, idx)
        assert_(np.allclose(B, A[:,idx[:k]], eps))
        assert_(np.allclose(B.dot(P), A, eps))

        # test SVD routines
        _debug_print("-----------------------------------------")
        _debug_print("SVD routines")
        _debug_print("-----------------------------------------")

        # fixed precision
        _debug_print("Calling iddp_svd / idzp_svd ...",)
        t0 = time.clock()
#.........这里部分代码省略.........

 def time_hilbert(self, size):
     sl.hilbert(size)

np.linalg.cond(A)
np.linalg.cond(A_delta)
b=np.matrix('1;1;1') #noise free right hand side


l=.01
k=100
ip.invert(A_delta,b,k,l)
np.linalg.pinv(A_delta)*b

#####################
#Example 3 Hilbert A
#####################

from scipy.linalg import hilbert
A=hilbert(3)
np.linalg.cond(A)
x=np.matrix('1;1;1')
b=A*x

#Now add noise to b of size ro
# dimension of A
r= np.matrix(A.shape)[0,0]
nb= np.random.normal(0, .50, r) #compute vector of normal variants mean=0,std=0.1
#nb= np.random.normal(0, 50, r) #compute vector of normal variants mean=0,std=50
#note that nb is 3 by 1, so we need to flip it
be=b+np.matrix(nb).transpose() #add noise

l=1
k=100