numpy

基础

from numpy import *

# 多项式的根：x^2 + 2x + 1
result = roots([1,2,1])

# 多元一次方程组
a = [[1,2],[2,-1]]
b = [4,3]
result = linalg.solve(a,b)

# 变形
x = arange(6)
x = x.reshape((2,3))

# 单位矩阵和零矩阵
eye(3)
zeros(9)
ones((2,6))

# 生成等间距样本数据
x = linspace(0,pi,10)

基础如下：

import numpy as np
"""
    备注：本示例中所有的matrix都使用的是array来体现
    实际上python中存在专门的matrix

    # 主要注意matrix和array的区别
    # matrix 使用更加接近matlab的使用方式
    简单示例如下：
    B = np.matrix('1,2;3,4')
    print(B*B)
    # 注意：若使用matrix，则*直接表示矩阵乘法
    # 而元素对应乘积则使用 multiply() 函数
"""
def vector_dot():
    a = np.array([1, 3, 5, 7, 9])
    b = np.array([1, 2, 3, 4, 5])
    # 对应元素相乘
    print(a * b)
    # 对应元素2次幂
    print(a ** 2)
    # dot product
    print(a.dot(b))
    print(type(a))
    print(a.dtype)
    print(a.shape)

def array_matrix_sum():
    d = np.arange(10)
    # 求和：针对数组
    print(d.sum())
    # 求最值
    print(np.min(d))    # 两种方式结果一致
    print(d.max())
    # 累加
    print(d.cumsum())

    x = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
    # 求和：针对矩阵
    print(x.sum())
    print(np.sum(x, axis=0))  # 列sum
    print(x.sum(axis=1))  # 行sum

def array_matrix_slice():
    # start,stop,num
    d = np.linspace(0, 90, 10)
    print(d)
    # start,stop,num,endpoint=False
    d = np.linspace(0, 100, 10, endpoint=False)
    print(d)
    # 切片
    indices = [1, 3, -1]
    print(d[indices])
    print(d[[1, 3, -1]])  # 外[]表示d的下标，内[]表示下标列表
    print(d[:5])
    # 按照条件 切片
    print(d[d >= 50])
    # 按照条件 查找元素下标
    index = np.where(d >= 50)
    print(index[0])
    print(d[index[0]])

    e = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12],
                  [13, 14, 15, 16]])
    print(e)
    # 矩阵切片方式
    print(e[1:3, 1:3])
    # 选取矩阵中某点
    print(e[3, 2])
    print(e[(3, 2)])
    # 花式索引 --- array[[x],[y]]
    # 坐标对应(1,1),(1,2),(2,1),(2,2)
    # 最终将回去的所有点排成一个一维数组
    print(e[[1, 1, 2, 2], [1, 2, 1, 2]])
    #
    print(e[np.arange(1, 3), 1:3])
    # 返回布尔值，（针对对应元素）
    # 该布尔值也可以用于索引切片
    print(e % 2 == 0)
    print(e[e % 2 == 0])

def array_stack():
    a = np.array([[1, 2, 3, 4], [5, 6, 7, 8]])
    b = np.array([0, 1, 0, -1])
    c = np.array([0, 1, 0, -1]).reshape(2, 2)
    print('-------------堆叠--------------')
    # 垂直堆叠
    r = np.vstack([a, b])
    print(r)
    # 水平堆叠
    r2 = np.hstack([a, c])
    print(r2)

    print('---------堆叠（方式2）----------')
    # 方式2（以垂直为例）
    b = np.array([[0, 1, 0, -1]])  # concatenate要求dimension一致
    r3 = np.concatenate((a, b), axis=0)
    # r3 = np.concatenate((a, c), axis=1)
    print(r3)

    print('---------堆叠（方式3）----------')
    # 方式3（以水平为例）
    # r4 = np.r_[a, b]  # 垂直堆叠
    r4 = np.c_[a, c]
    print(r4)

    print('-------------扩展--------------')
    x = np.array([1, 2, 3])
    # repeat以元素为单位扩展
    ret = np.repeat(x, 3)
    print(ret)
    # repeat以数组为单位扩展
    ret2 = np.tile(x, 3)
    print(ret2)
    ret3 = np.r_[np.repeat(x, 3), np.tile(x, 3)]
    print(ret3)

def create_matrix():
    a = np.zeros((3, 3))
    print(a)
    b = np.ones((3, 2))
    print(b)
    # 用9填充3*3的矩阵
    c = np.full((3, 3), 9)
    print(c)
    # 用6填充类似b的矩阵
    d = np.full_like(b, 6)
    print(d)
    e = np.eye(3)
    print(e)
    r = np.random.random((2, 2))
    print(r)

def dot_product():
    # 针对 1 dimension
    # 此时dot作用类似于inner(内积)
    # 作用式：两向量数量积
    v = np.array([9, 10])
    w = np.array([11, 12])
    print(v.dot(w))         # Two different ways to express
    print(np.dot(v, w))     # ret = 219

    # 针对 2 dimension
    # 作用是：矩阵乘法
    x = np.array([[1, 2], [3, 4]])
    y = np.array([[5, 6], [7, 8]])
    print(np.dot(x, y))  # dot for 2d -> Matrix / vector product
    # ret = [[19,22],[43,50]]

    # 针对 矩阵
    # 作用是：矩阵乘法
    x = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
    print(x)
    y = np.eye(3)
    y[(1, 0)] = 2   # 用矩阵描述 矩阵消元法（高斯消元法）
    print(y)
    print(np.dot(y, x))

def matrix_transpose():
    x = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
    print(x.T)
    # @attention 1d   .T无效
    y = np.array([1, 2, 3])
    print(y.T)

    # 翻转行&列的方式
    print(x[::-1])          # 翻转行
    print(x[:, ::-1])       # 翻转列

def broadcasting(mode=2):
    """
        默认表示经典用法。其他表示常规用法
    """
    x = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12]])
    v = np.array([1, 0, -1])
    if mode == 0:
        """
            当数据大时，效率低下
        """
        y = np.empty_like(x)
        for i in range(4):
            y[i, :] = x[i, :] + v
    elif mode == 1:
        vv = np.tile(v, (4, 1))
        y = x + vv
    elif mode == 2:
        """
            广播
        """
        y = x + v
    else:
        y = None
    print(y)

    """
        广播应用2
    """
    a = np.array([1, 2, 3])
    b = np.array([4, 5])
    c = a.reshape((3, 1))
    print(c)
    print(c * b)

def linalg_demo():
    # 资料链接
    # https://www.numpy.org.cn/article/basics/numpy_matrices_vectors.html

    print('----- 解线性方程组 -----')
    A = np.array([[2, 1, -2], [3, 0, 1], [1, 1, -1]])
    b = np.transpose(np.array([[-3, 5, -2]]))
    x = np.linalg.solve(A, b)
    print(x)
    # 查看一些结果，看是否正确
    print(A.dot(x))

    print('----- 求矩阵的逆 -----')
    C = np.linalg.inv(A)
    # 使用round的目的是，让特别小的数，例如1.11e-16显示为0
    print(C)
    print(np.round(C, 6))
    ret = C.dot(A)
    print(np.round(ret, 6))

    print('----- 计算行列式 -----')
    # 方式1
    # D = np.array([[3, 4], [5, 6]])
    # E = np.asmatrix(D)
    # 方式2
    E = np.mat('3,4; 5,6')
    print(E)
    ret = np.linalg.det(E)
    print(np.round(ret, 6))

def float_random_demo():
    a = np.random.random((10, 10))
    np.set_printoptions(precision=3)  # 显示精度3位小数
    np.set_printoptions(threshold=3)  # 最多显示6项，其他项使用省略号
    # np.set_printoptions(threshold=np.nan)   # 显示所有项目
    print(a)

def uniformization():
    value = np.random.randint(5, 10, size=16).reshape(4, 4)
    print(value)
    v_max, v_min = value.max(), value.min()
    # 归一化处理
    uniform_value = (value - v_min) / (v_max - v_min)
    print(uniform_value)

if __name__ == "__main__":
    # vector_dot()
    # array_matrix_sum()
    # array_matrix_slice()
    # array_stack()
    # create_matrix()
    # dot_product()
    # matrix_transpose()
    # broadcasting(0)
    # linalg_demo()
    # float_random_demo()
    uniformization()

点乘

from scipy.spatial.distance import pdist
import numpy as np

def get_module(x):
    """
        x为需要求长度的向量
        欧拉距离/欧式距离(Euclidean Distance)
    """
    return np.sqrt(np.sum(np.square(x)))

def get_angle(x, y):
    """
        一维数据，dot点乘，对数学而言即为求取其数量积
    """
    return np.arccos(np.dot(x, y)/(get_module(x)*get_module(y)))

def get_mudule_len_by_dot():
    print('----- 通过数量积求模长 ----- ')
    # --------------------------------
    # 纯numpy(数学)方式
    OA = np.array([3, 0])
    OB = np.array([0, 4])
    AB = OB - OA  # 指向被减数
    # 长度
    d = get_module(AB)
    print('mlen = ', d)

    # -------------------------------
    # numpy + scipy方式
    xy = np.vstack([OA, OB])
    print(xy)
    dist = pdist(xy)
    print('mlen = ', dist[0])

def get_angle_by_cross():
    """
        已知M(1,1,1),A(2,2,1),B(2,1,2)
        求角度AMB
    """
    print('----- 通过向量积求角度 ----- ')
    M = np.array([1, 1, 1])
    A = np.array([2, 2, 1])
    B = np.array([2, 1, 2])
    MA = A - M
    MB = B - M
    angle = get_angle(MA, MB)
    print(angle)
    if angle - np.pi/3 < 0.0001:
        print('数量积结果正确')

def get_triangle_area_by_cross():
    """
        已知三角形三点坐标，A，B，C坐标点。
        求三角形面积
    """
    print('----- 通过向量积求三角形面积 ----- ')
    A = np.array([1, 2, 3])
    B = np.array([3, 4, 5])
    C = np.array([2, 4, 7])
    AB = B - A
    AC = C - A
    # 面积
    AB_AC = 1/2 * np.cross(AB, AC)
    S = get_module(AB_AC)
    print(S)
    if S - np.sqrt(14) < 0.0001:
        print('向量积结果正确')

if __name__ == "__main__":
    get_mudule_len_by_dot()
    get_angle_by_cross()
    get_triangle_area_by_cross()

叉乘

import numpy as np

def coor_trans(point, theta):
	"""
		coordinate transformation (坐标转换)

		theta方向：以顺时针旋转为正
	"""
	point = np.transpose(point)
	k = np.array([[np.cos(theta), np.sin(theta)],
				  [-np.sin(theta), np.cos(theta)]])
	print(point)
	# return np.dot(k, point)
	return np.round(np.dot(k, point),6)


if __name__ == "__main__":
	point = np.array([[1,0],[1,2],[2,0],[1,0]])
	print(point)
	point_trans = coor_trans(point, np.pi / 2)
	print(point_trans)

矩阵变换

import numpy as np
import scipy.linalg

def matrix_transform():
    """
        目的是通过矩阵乘法实现变换 ------ 旋转，翻转
        **直接**进行切片操作的此处就不介绍了
    """
    # A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
    A = np.arange(1, 10).reshape(3, 3)
    print('----- 原始矩阵 -----\n', A)

    B = np.eye(3)[[2, 1, 0]]
    # C = np.transpose(A)
    C = A.T

    print('----- 上下翻转 ----')
    A1 = B.dot(A)       # 左乘
    print(A1)
    print('----- 左右翻转 ----')
    A2 = A.dot(B)       # 右乘
    print(A2)

    print('----- 旋转180° ----')
    print(np.dot(A1, B))    # rotate180 = filpud + filplr
    print(np.dot(B, A2))    # rotate180 = filplr + filpud

    print('----- 旋转90°(以顺时针为正) ----')
    print(np.dot(C, B))     # rotate90 = transpose + filpud
    print('----- 旋转-90°(以顺时针为正) ----')
    print(np.dot(B, C))     # rotate-90 = transpose + filplr


def matrix_decomposition_lu():
    """
        https://baike.baidu.com/item/lu%E5%88%86%E8%A7%A3/764245?fr=aladdin
        scipy采用'列主元三角分解': A=PLU
        手动计算(高斯消元法): A=LU
    """
    A = np.matrix('2,3; 5,4')
    b = np.matrix('4,3').T

    print('----- LU分解 (返PLU参数)-----')
    p, l, u = scipy.linalg.lu(A)
    print(p)
    print(l)
    print(u)

    print('----- Ax=b (A=PLU) 求解x -----')
    y = np.linalg.solve(np.dot(p, l), b)
    x = np.linalg.solve(u, y)
    print(x)

    print('----- 直接解方程 -----')
    x = np.linalg.solve(A, b)
    print(x)

    """
    # 这个L相当于是上述的PL，形式很怪，此处不采用
    print('----- LU分解 (返LU参数)-----')
    l, u = scipy.linalg.lu(A, permute_l=True)
    print(l)
    print(u)
    """


if __name__ == "__main__":
    # matrix_transform()
    matrix_decomposition_lu()

卷积

import numpy as np

'''
	py卷积计算
'''
a = np.array([1,2])
b = np.array([3,4])
print(a)
print(b)
ret = np.convolve(a,b)
print(ret)

'''
	系数向量内积
'''
c = np.array([[2,1,0,],[0,2,1],[0,0,2]])
d = np.array([0,3,4])
print(c)
print(d)
ret = c.dot(d)
# print(ret)


def lsl(lst,k,cnt):
	temp = np.concatenate((lst[-k:],lst[:-k]),axis=0)
	z = np.zeros(cnt)
	z[k:] = temp[k:]
	return z

def conv(a,b):
	'''
		将卷积转换为矩阵运算
	'''
	l_a = len(a)
	l_b = len(b)
	item_cnt = len(a) + len(b) - 1
	data = np.zeros((item_cnt,item_cnt))
	data2 = np.zeros(item_cnt)
	data2[-l_b:] = b[-l_b:]
	a = a[::-1]
	for i in range(item_cnt):
		data[i][:l_a] = a
		data[i][:] = lsl(data[i],i,item_cnt)
	return np.dot(data,data2.transpose())


ret = conv(a,b)
print(ret)

最小二乘法

# -*- coding: utf-8 -*-
"""
Created on Sun Nov 17 14:12:55 2019

@author: Administrator
"""

import matplotlib.pyplot as plt
import numpy as np
x = np.array([25,27,31,33,35])
y = np.array([110,115,155,160,180])
# result = (7.2,-73)
X = np.mat([[25,1],[27,1],[31,1],[33,1],[35,1]])
Y = np.mat(y).T
"""
    最小二乘法的矩阵形式
"""
X_T = X.transpose()
A = np.dot(np.dot(np.linalg.inv(np.dot(X_T,X)),X_T),Y)

a = A[0,0]
b = A[1,0]
y2 = a*x + b
plt.plot(x,y2)
plt.scatter(x,y)
plt.show()

梯度下降法

# -*- coding: utf-8 -*-
"""
Created on Mon Nov 18 20:46:20 2019

@author: Administrator
"""

import matplotlib.pyplot as plt
import numpy as np

def cost_function(X,y,theta):
    m = len(y)
    y_mat = y.reshape(m,1)
    J = 1/(2*m) * np.sum(np.power((X@theta - y_mat),2))
    return J

def gradient_descent(X,y,theta,alpha,num_iters):
    J_history = np.zeros((num_iters,1))
    m = len(y)
    y_mat = y.reshape(m,1)
    for n_iter in range(num_iters):
        theta = theta - alpha * 1 / m * X.T @ (X @ theta - y_mat)
        J_history[n_iter, 0] = cost_function(X,y,theta)
    return J_history, theta

"""
the result of least square = 
[[  7.20930233]
 [-73.72093023]]
"""
x = np.array([25,27,31,33,35])
y = np.array([110, 115, 155, 160, 180])
m = len(x)
X = np.hstack((x.reshape(m, 1), np.ones((m, 1))))

# initial value = (0,0)
theta = np.array([0, 0]).reshape(2, 1)
# if alpha >= 0.0022 ---> overflow (@.@)
alpha = 0.0021

J,theta = gradient_descent(X,y,theta,alpha,400000)
print(theta)