高级编程技术第十一周作业

最新推荐文章于 2018-06-13 10:51:51 发布

Dexan_N

最新推荐文章于 2018-06-13 10:51:51 发布

阅读量202

点赞数

本文链接：https://blog.csdn.net/Dexan_N/article/details/80411928

版权

本周作业是numpy库的使用。这个库在python环境下提供了许多数学相关的函数及类，以期望在python环境下实现matlab的计算。不过既然是python，效率也不要太强求了。

Exercise 9.1: Matrix operations

Calculate A + A, AA>;A>A and AB. Write a function that computes A(B - λI) for any .

当然在题目给定的条件下，n=200、m=500实在过于巨大（文本量而言），导致基本没法将结果记录下来，因此下面统一采用n=4、m=10。

import random
import numpy
import numpy.matlib
from scipy.linalg import toeplitz


def f(a, b, m, lda):
    i = numpy.matlib.identity(m) * lda
    return numpy.dot(a, b - i)


n = 4
m = 10
a = numpy.matlib.empty((n, m))
for i in range(0, n * m):
    a.flat[i] = random.gauss(0, 1)
b = toeplitz(range(1, m + 1))
print("A + A:")
print(a + a)
print("A * At:")
print(numpy.dot(a, numpy.transpose(a)))
print("At * A:")
print(numpy.dot(numpy.transpose(a), a))
print("A * B:")
print(numpy.dot(a, b))
lda = input("Please input the lambda:")
print("A * (B - lambda * I):")
print(f(a, b, m, float(lda)))

结果如下：

A + A:
[[-0.84938612 -1.1561019   2.75438702 -1.71456328  0.04059118  0.1534713
  -1.66728747 -0.01302837 -1.34270901 -1.87653223]
 [-0.45028936 -0.6618186  -2.81645558 -2.64948963  1.12849407  2.64990658
  -5.05316526 -3.88338383  1.34767272  1.62100678]
 [-1.78059658  0.63376315 -3.33167737 -0.49458535  0.28742315  2.85498364
   1.06723189 -0.31354899 -2.30749876  3.98076355]
 [-1.27076995  1.89063576  1.61928496  2.97298097 -0.82483351  2.33784495
  -1.87384472 -0.39021517 -2.93498161 -0.18986124]]
A * At:
[[ 5.17846563  0.50236731 -3.31155399  1.50203002]
 [ 0.50236731 17.23687984  4.53347157 -0.28281946]
 [-3.31155399  4.53347157 11.38968454  1.79306182]
 [ 1.50203002 -0.28281946  1.79306182  8.77741375]]
At * A:
[[ 1.42739943e+00 -5.62763205e-01  7.00829136e-01 -6.19886627e-02
  -1.55937670e-03 -2.34450458e+00  1.04311837e+00  7.03472626e-01
   2.09301164e+00 -1.49572142e+00]
 [-5.62763205e-01  1.43768369e+00 -9.25964160e-02  2.26076639e+00
  -5.42771856e-01  1.07455269e+00  6.01361844e-01  4.12171806e-01
  -1.58774883e+00  8.15138205e-01]
 [ 7.00829136e-01 -9.25964160e-02  7.31030694e+00  2.30037529e+00
  -1.33994783e+00 -3.19171801e+00  7.62424904e-01  2.82856684e+00
  -1.13968270e+00 -5.82606223e+00]
 [-6.19886627e-02  2.26076639e+00  2.30037529e+00  4.76068826e+00
  -1.41347477e+00 -4.36425203e-01  2.53705922e+00  2.32657440e+00
  -2.21321869e+00 -9.02672145e-01]
 [-1.55937670e-03 -5.42771856e-01 -1.33994783e+00 -1.41347477e+00
   5.09527223e-01  4.72222251e-01 -9.79446783e-01 -1.03779079e+00
   8.05995285e-01  7.63473353e-01]
 [-2.34450458e+00  1.07455269e+00 -3.19171801e+00 -4.36425203e-01
   4.72222251e-01  5.16550222e+00 -3.74503137e+00 -3.02501092e+00
  -2.52106590e+00  3.73216784e+00]
 [ 1.04311837e+00  6.01361844e-01  7.62424904e-01  2.53705922e+00
  -9.79446783e-01 -3.74503137e+00  8.24115115e+00  5.01041887e+00
  -3.83566892e-01 -1.14582047e-01]
 [ 7.03472626e-01  4.12171806e-01  2.82856684e+00  2.32657440e+00
  -1.03779079e+00 -3.02501092e+00  5.01041887e+00  3.83285514e+00
  -8.36812224e-01 -1.86115525e+00]
 [ 2.09301164e+00 -1.58774883e+00 -1.13968270e+00 -2.21321869e+00
   8.05995285e-01 -2.52106590e+00 -3.83566892e-01 -8.36812224e-01
   4.38943921e+00 -9.81036080e-01]
 [-1.49572142e+00  8.15138205e-01 -5.82606223e+00 -9.02672145e-01
   7.63473353e-01  3.73216784e+00 -1.14582047e-01 -1.86115525e+00
  -9.81036080e-01  5.50789048e+00]]
A * B:
[[-21.62892049 -19.64272717 -18.81263575 -15.2281573  -13.35824214
  -11.44773579  -9.38375815  -8.98706797  -8.60340616  -9.56245337]
 [-18.68972355 -14.75625186 -11.48459877 -11.02940126 -13.22369337
  -14.28949142 -12.70538288 -16.1744396  -23.52688015 -29.53164798]
 [ 15.04147591  12.96275016  11.51778756   6.74114759   1.46992227
   -3.5138799   -5.64269843  -6.70428507  -8.0794207  -11.76205508]
 [ -7.69444942  -9.63333959  -9.68159399  -8.11056344  -3.56655192
    0.15262609   6.20964905  10.3928273   14.18579037  15.04377183]]
Please input the lambda:0.5
A * (B - lambda * I):
[[-21.41657396 -19.35370169 -19.5012325  -14.79951648 -13.36838993
  -11.48610362  -8.96693628  -8.98381087  -8.26772891  -9.09332031]
 [-18.57715121 -14.59079721 -10.78048488 -10.36702885 -13.50581689
  -14.95196806 -11.44209156 -15.20359364 -23.86379833 -29.93689968]
 [ 15.48662505  12.80430938  12.35070691   6.86479393   1.39806648
   -4.22762581  -5.90950641  -6.62589783  -7.50254601 -12.75724597]
 [ -7.37675693 -10.10599853 -10.08641523  -8.85380869  -3.36034355
   -0.43183515   6.67811023  10.49038109  14.91953577  15.09123714]]

Exercise 9.2: Solving a linear system

Generate a vector b with m entries and solve Bx = b.

可以通过求逆来求解x，当然numpy直接有求解函数，就直接用吧。。。

import numpy
import numpy.matlib
import numpy.linalg
from scipy.linalg import toeplitz

n = 4
m = 10
B = toeplitz(range(1, m + 1))
b = numpy.matlib.rand((m, 1))
x = numpy.linalg.solve(B, b)
print("B =")
print(B)
print("b = ")
print(b)
print("x = ")
print(x)
print("Bx = ")
print(numpy.dot(B, x))

结果如下：

B =
[[ 1  2  3  4  5  6  7  8  9 10]
 [ 2  1  2  3  4  5  6  7  8  9]
 [ 3  2  1  2  3  4  5  6  7  8]
 [ 4  3  2  1  2  3  4  5  6  7]
 [ 5  4  3  2  1  2  3  4  5  6]
 [ 6  5  4  3  2  1  2  3  4  5]
 [ 7  6  5  4  3  2  1  2  3  4]
 [ 8  7  6  5  4  3  2  1  2  3]
 [ 9  8  7  6  5  4  3  2  1  2]
 [10  9  8  7  6  5  4  3  2  1]]
b =
[[0.19485856]
 [0.06491564]
 [0.45953234]
 [0.35550435]
 [0.80744295]
 [0.14914871]
 [0.43569777]
 [0.95209637]
 [0.30585147]
 [0.86331063]]
x =
[[-0.01687286]
 [ 0.26227981]
 [-0.24932235]
 [ 0.2779833 ]
 [-0.55511642]
 [ 0.47242165]
 [ 0.11492477]
 [-0.58132175]
 [ 0.60185203]
 [-0.23063098]]
Bx =
[[0.19485856]
 [0.06491564]
 [0.45953234]
 [0.35550435]
 [0.80744295]
 [0.14914871]
 [0.43569777]
 [0.95209637]
 [0.30585147]
 [0.86331063]]

Exercise 9.3: Norms

Compute the Frobenius norm of A: ||A||F and the innity norm of B: ||B||1. Also nd the largest and

smallest singular values of B.

仍旧使用自带的函数即可

import random
import numpy
import numpy.matlib
import numpy.linalg
from scipy.linalg import toeplitz

n = 4
m = 10
a = numpy.matlib.empty((n, m))
for i in range(0, n * m):
    a.flat[i] = random.gauss(0, 1)
b = toeplitz(range(1, m + 1))
print("A:")
print(a)
print("B:")
print(b)
print("Frobenius norm of A:")
print(numpy.linalg.norm(a, 'fro'))
print("innity norm of B:")
print(numpy.linalg.norm(b, numpy.inf))
print("largest singular values of B:")
print(numpy.linalg.norm(b, 2))
print("smallest singular values of B:")
print(numpy.linalg.norm(b, -2))

结果如下：

A:
[[ 0.00270799  1.55158354  1.4929957   1.87874988  0.06912358 -1.33506174
  -0.85171565 -0.59300652  1.67936931 -0.19606906]
 [ 0.65661568 -0.22044142 -1.63375831 -0.08908255  0.98427599  1.15290236
  -0.01545519  0.25548508 -1.82891623 -1.51978535]
 [ 0.1118962  -0.50953158 -0.42362632 -1.09725616 -0.18091825 -0.5525389
   0.20537803  0.3827813  -0.73509675 -2.3665386 ]
 [-0.5988766  -0.82232487  1.45208276  1.97098965 -0.71494838  0.03329491
   1.8172782   1.99291261  1.04028339 -0.0037711 ]]
B:
[[ 1  2  3  4  5  6  7  8  9 10]
 [ 2  1  2  3  4  5  6  7  8  9]
 [ 3  2  1  2  3  4  5  6  7  8]
 [ 4  3  2  1  2  3  4  5  6  7]
 [ 5  4  3  2  1  2  3  4  5  6]
 [ 6  5  4  3  2  1  2  3  4  5]
 [ 7  6  5  4  3  2  1  2  3  4]
 [ 8  7  6  5  4  3  2  1  2  3]
 [ 9  8  7  6  5  4  3  2  1  2]
 [10  9  8  7  6  5  4  3  2  1]]
Frobenius norm of A:
7.020270644384
innity norm of B:
55.0
largest singular values of B:
44.076292527938
smallest singular values of B:
0.5125428154684579

Exercise 9.4: Power iteration

Generate a matrix Z, n x n, with Gaussian entries, and use the power iteration to nd the largest
eigenvalue and corresponding eigenvector of Z. How many iterations are needed till convergence?

Optional: use the time.clock() method to compare computation time when varying n.

使用幂迭代法求特征值需要有一个比较好的初值，否则很容易会算不出来。然而由于Z本身是高斯分布的（也就是说它本身可能不收敛），所以特征值干脆也用随机，让RNG来解决一切。

import random
import numpy
import numpy.matlib
import numpy.linalg

n = 10
z = numpy.matlib.empty((n, n))
for i in range(0, n * n):
    z.flat[i] = random.gauss(0, 1)
b = numpy.random.rand(z.shape[1], 1)
ck = True
times = 0
while(ck):
    bak = b
    ck = True
    times = times + 1
    bk = numpy.dot(z, b)
    bk_norm = numpy.linalg.norm(bk)
    b = bk / bk_norm
    if(numpy.linalg.norm(b - bak, numpy.inf) < 1 / 2 * 10 * pow(10, -5)):
        ck = False
print("eigenvector:")
print(b)
print("eigenvalue:")
print(numpy.linalg.norm(numpy.dot(z, b) / b))
print("iteration times:")
print(times)

某次成功的结果如下：

eigenvector:
[[ 0.42296718]
 [ 0.08746957]
 [ 0.10301223]
 [-0.16572223]
 [ 0.6286509 ]
 [ 0.04778553]
 [-0.05736135]
 [ 0.3854316 ]
 [-0.37260999]
 [-0.29529816]]
eigenvalue:
13.314243915992803
iteration times:
65

Exercise 9.5: Singular values

Generate an nn matrix, denoted by C, where each entry is 1 with probability p and 0 otherwise. Use
the linear algebra library of Scipy to compute the singular values of C. What can you say about the

relationship between n, p and the largest singular value?

可以用随机分布来生成C，然后使用scipy求无穷范数（通常为最大特征值）。

import random
import numpy
import numpy.matlib
import numpy.linalg
from scipy.linalg import eig

n = 10
p = 0.6
a = numpy.matlib.zeros((n, n))
for i in range(0, n * n):
    if(random.random() < p):
        a.flat[i] = 1
w, v = eig(a)
print("all singular value:")
print(w)
print("the largest singular value:")
print(numpy.linalg.norm(w, numpy.inf))

某次的运行结果：

all singular value:
[ 6.3127503 +0.j         -0.92315776+0.22821323j -0.92315776-0.22821323j
 -0.81560421+0.57041439j -0.81560421-0.57041439j  0.38047453+1.07777556j
  0.38047453-1.07777556j  1.79216704+0.j          0.36292569+0.j
  1.24873183+0.j        ]
the largest singular value:
6.312750304338822

多次运行可以看到，最大特征值的均值大致等同于np。

Exercise 9.6: Nearest neighbor

Write a function that takes a value z and an array A and nds the element in A that is closest to z. The
function should return the closest value, not index.
Hint: Use the built-in functionality of Numpy rather than writing code to nd this value manually. In

particular, use brackets and argmin.

argmin可以求出整个矩阵/数组中最小元的一维迭代坐标，然而若是来求和某个数值最近的元则可以通过让整个矩阵/数组和这个数值相减并取绝对值来解决。

import numpy


def f(z, A, n):
    ak = A - z
    for i in range(0, n):
        ak.flat[i] = abs(ak.flat[i])
    return a.flat[numpy.argmin(ak)]


n = 10
a = numpy.random.rand(n) * 100
z = 50
print("A:")
print(a)
print("z = ")
print(z)
print("the cloest value:")
print(f(z, a, n))

结果如下：

A:
[96.54910492 23.62329268 85.24039685 46.90811694 96.08257854 29.16258579
  5.72084394 94.55313671 15.41269397 26.13132433]
z =
50
the cloest value:
46.908116935194165