向量
- 线性代数是从研究一个数拓展到一组数
- 一组数的基本表示方法——向量(Vector)
- 向量是线性代数研究的基本元素
- 一组数的作用：最基本的出发点：表示方向
在二维空间中，不同的方向表示不同的终点，比如同样走5000米
它的终点是(4,3)
而它的终点是(4.6,1)
在三维空间中也是如此
在线性代数的世界里，起始点不重要
在这个图中，从(-1,-1)到(3,2)和从(0,0)到(4,3)是一样的。它们只是坐标系不同而已。
为了研究方便，我们定义向量都从原点起始，但是顺序是重要的，很显然(4,3)和(3,4)是不同的。
向量是一组有序的数。
- 如果只是表示方向，最多三个维度就够了，但为了扩展向量所能表达的范围，更加抽象的：n维向量。比如一个房子的维度
面积 卧室 卫生间 最近地铁站(km) 价格(万)
120 3 2 2 666
所以这是一个五维的向量(120,3,2,2,666),此时向量就是一组数，这组数的含义由使用者定义。
- 两个视角看似不同，但可以互相转换，比如房子的这个向量，我们可以想象它是一个五维空间的一个终点。
- 而一个方向就是一个点
- 空间中的一个点，可以看作从原点指向这个点的一个方向
- 在学习初始，使用方向的视角，更直观，更形象。
- 更关键的是：这两个视角，都不是简单的"一组数"
- 一个是一个有向线段
- 一个是空间中的点
更严格一些的定义
- 和向量相对应，一个数字，称为标量
- 代数，用符号代表数。和标量相区别，向量的符号画箭头:
- 个别情况下，尤其是几何学中，我们会考虑向量的起始点
行向量和列向量 (3,4) 
通常教材，论文，提到向量，都指列向量。由于横版印刷，使用符号：
，这个T是一个运算符号，把这个向量进行转置的操作。
向量的类，因为向量就是一组数，所以我们用数组来表示。
接口
public interface Attribute<T extends Number> { /** * 向量的长度 * @return */ int len(); /** * 获取向量的元素 * @param index 索引 * @return */ T get(int index); }
实现类
public class Vector<T extends Number> implements Attribute<T> { private T[] values; public Vector(T[] values) { this.values = values; } @Override public String toString() { return "Vector{" + "values=" + Arrays.toString(values) + '}'; } @Override public int len() { return values.length; } @Override public T get(int index) { if (index >= len() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } return values[index]; } public static void main(String[] args) { Attribute<Integer> vector = new Vector<>(new Integer[]{5,2}); System.out.println(vector.toString()); System.out.println(vector.len()); System.out.println(vector.get(0)); } }
运行结果
Vector{values=[5, 2]}
2
5
Python代码
class Vector: def __init__(self,lst): self._values = lst def __getitem__(self, index): # 取向量的第index个元素 return self._values[index] def __len__(self): # 返回向量的长度 return len(self._values) def __repr__(self): return "Vector({})".format(self._values) def __str__(self): return "({})".format(", ".join(str(e) for e in self._values))
from playLA.Vector import Vector if __name__ == "__main__": vec = Vector([5,2]) print(vec) print(len(vec)) print("vec[0] = {}, vec[1] = {}".format(vec[0],vec[1]))
运行结果
(5, 2)
2
vec[0] = 5, vec[1] = 2
向量的基本运算
- 向量加法
向量的加法类似于我们中学物理中的矢量的平行四边形法则，它相当于先向(5,2)的点走一段距离，再从终点向(2,5)的方向再走一段距离。最终得到的地点就是(7,7)。
所以一个二维向量的加法就是
三维向量同理
N维向量同理
- 数量乘法
向量的乘法其实就是在向量的同一个方向上再移动乘数倍的距离，比如这里是2倍，就会移动2倍距离，得到终点为(10,4)
所以向量乘法就是
N维向量同理
在编码中，我们需要扩展Attribute接口
public interface Attribute<T extends Number> { /** * 向量的长度 * @return */ int len(); /** * 获取向量的元素 * @param index 索引 * @return */ T get(int index); /** * 向量相加 * @param another * @return */ Attribute<T> add(Attribute<T> another,Addr<T> addr); /** * 向量乘以一个数 * @param k * @param addr * @return */ Attribute<T> mul(T k,Addr<T> addr); }
新增一个Addr的接口，不过该接口没有实现类，而是在实际编码中去编写匿名实现类
public interface Addr<T extends Number> { T zero(); T add(T lhs,T rhs); T mul(T k,T hs); }
最后是向量的实现类
public class Vector<T extends Number> implements Attribute<T> { private T[] values; public Vector(T[] values) { this.values = values; } @Override public String toString() { return "Vector{" + "values=" + Arrays.toString(values) + '}'; } @Override public int len() { return values.length; } @Override public T get(int index) { if (index >= len() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } return values[index]; } @Override @SuppressWarnings("unchecked") public Attribute<T> add(Attribute<T> another,Addr<T> addr) { if (this.len() != another.len()) { throw new IllegalArgumentException("加法错误，两个向量的维度必须相等"); } T[] newValues = (T[])new Number[this.len()]; for (int i = 0; i < this.len(); i++) { List<T> list = new ArrayList<>(); list.add(this.get(i)); list.add(another.get(i)); newValues[i] = sumValue(list, addr); } this.values = newValues; return this; } @Override @SuppressWarnings("unchecked") public Attribute<T> mul(T k, Addr<T> addr) { T[] newValues = (T[])new Number[this.len()]; for (int i = 0; i < this.len(); i++) { newValues[i] = addr.mul(k,this.get(i)); } this.values = newValues; return this; } private T sumValue(List<T> list, Addr<T> addr) { T total = addr.zero(); for (T n : list) { total = addr.add(total,n); } return total; } public static void main(String[] args) { Attribute<Integer> vector = new Vector<>(new Integer[]{5,2}); System.out.println(vector.toString()); System.out.println(vector.len()); System.out.println(vector.get(0)); Attribute<Integer> vector1 = new Vector<>(new Integer[]{6,3}); Addr<Integer> addr = new Addr<Integer>() { @Override public Integer zero() { return 0; } @Override public Integer add(Integer lhs, Integer rhs) { return lhs + rhs; } @Override public Integer mul(Integer k, Integer hs) { return k * hs; } }; vector.add(vector1, addr); System.out.println(vector); vector.mul(3,addr); System.out.println(vector); } }
运行结果
Vector{values=[5, 2]}
2
5
Vector{values=[11, 5]}
Vector{values=[33, 15]}
Python代码
class Vector: def __init__(self, lst): self._values = lst def __add__(self, another): # 向量加法，返回结果向量 assert len(self) == len(another),\ "加法错误，两个向量的维度必须相等" return Vector([a + b for a, b in zip(self, another)]) def __mul__(self, k): # 返回数量乘法的结果向量 return Vector([k * e for e in self]) def __iter__(self): # 返回向量的迭代器 return self._values.__iter__() def __getitem__(self, index): # 取向量的第index个元素 return self._values[index] def __len__(self): # 返回向量的长度 return len(self._values) def __repr__(self): return "Vector({})".format(self._values) def __str__(self): return "({})".format(", ".join(str(e) for e in self._values))
from playLA.Vector import Vector if __name__ == "__main__": vec = Vector([5, 2]) print(vec) print(len(vec)) print("vec[0] = {}, vec[1] = {}".format(vec[0], vec[1])) vec2 = Vector([6, 3]) print("{} + {} = {}".format(vec, vec2, vec + vec2)) print("{} * {} = {}".format(vec, 3, vec * 3))
运行结果
(5, 2)
2
vec[0] = 5, vec[1] = 2
(5, 2) + (6, 3) = (11, 5)
(5, 2) * 3 = (15, 6)
向量运算的基本性质
加法交换律 
加法结合律 
乘法分配律 
乘法交换律 
零向量
对于任意一个向量
,都存在一个向量O,满足 + O = ，根据向量的加法，我们可以推导出
由此我们可以得到零向量的各个维度的值都是0.并且零向量O没有箭头。
同此，我们可以推导，对于任意一个向量，都存在一个向量
，满足+ = O。上述唯一。
我们给出一个创建零向量的静态方法
public class Vector<T extends Number> implements Attribute<T> { private T[] values; public Vector(T[] values) { this.values = values; } @SuppressWarnings("unchecked") public static <T extends Number> Attribute<T> zero(int dim,Addr<T> addr) { T[] values = (T[])new Number[dim]; Arrays.fill(values,addr.zero()); return new Vector<>(values); } @Override public String toString() { return "Vector{" + "values=" + Arrays.toString(values) + '}'; } @Override public int len() { return values.length; } @Override public T get(int index) { if (index >= len() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } return values[index]; } @Override @SuppressWarnings("unchecked") public Attribute<T> add(Attribute<T> another,Addr<T> addr) { if (this.len() != another.len()) { throw new IllegalArgumentException("加法错误，两个向量的维度必须相等"); } T[] newValues = (T[])new Number[this.len()]; for (int i = 0; i < this.len(); i++) { List<T> list = new ArrayList<>(); list.add(this.get(i)); list.add(another.get(i)); newValues[i] = sumValue(list, addr); } this.values = newValues; return this; } @Override @SuppressWarnings("unchecked") public Attribute<T> mul(T k, Addr<T> addr) { T[] newValues = (T[])new Number[this.len()]; for (int i = 0; i < this.len(); i++) { newValues[i] = addr.mul(k,this.get(i)); } this.values = newValues; return this; } private T sumValue(List<T> list, Addr<T> addr) { T total = addr.zero(); for (T n : list) { total = addr.add(total,n); } return total; } public static void main(String[] args) { Attribute<Integer> vector = new Vector<>(new Integer[]{5,2}); System.out.println(vector.toString()); System.out.println(vector.len()); System.out.println(vector.get(0)); Attribute<Integer> vector1 = new Vector<>(new Integer[]{6,3}); Addr<Integer> addr = new Addr<Integer>() { @Override public Integer zero() { return 0; } @Override public Integer add(Integer lhs, Integer rhs) { return lhs + rhs; } @Override public Integer mul(Integer k, Integer hs) { return k * hs; } }; vector.add(vector1, addr); System.out.println(vector); vector.mul(3,addr); System.out.println(vector); System.out.println(Vector.zero(4,addr)); } }
运行结果
Vector{values=[5, 2]}
2
5
Vector{values=[11, 5]}
Vector{values=[33, 15]}
Vector{values=[0, 0, 0, 0]}
Python代码
class Vector: def __init__(self, lst): self._values = lst @classmethod def zero(cls,dim): # 返回一个dim维的零向量 return cls([0] * dim) def __add__(self, another): # 向量加法，返回结果向量 assert len(self) == len(another),\ "加法错误，两个向量的维度必须相等" return Vector([a + b for a, b in zip(self, another)]) def __mul__(self, k): # 返回数量乘法的结果向量 return Vector([k * e for e in self]) def __iter__(self): # 返回向量的迭代器 return self._values.__iter__() def __getitem__(self, index): # 取向量的第index个元素 return self._values[index] def __len__(self): # 返回向量的长度 return len(self._values) def __repr__(self): return "Vector({})".format(self._values) def __str__(self): return "({})".format(", ".join(str(e) for e in self._values))
from playLA.Vector import Vector if __name__ == "__main__": vec = Vector([5, 2]) print(vec) print(len(vec)) print("vec[0] = {}, vec[1] = {}".format(vec[0], vec[1])) vec2 = Vector([6, 3]) print("{} + {} = {}".format(vec, vec2, vec + vec2)) print("{} * {} = {}".format(vec, 3, vec * 3)) zero2 = Vector.zero(4) print(zero2)
运行结果
(5, 2)
2
vec[0] = 5, vec[1] = 2
(5, 2) + (6, 3) = (11, 5)
(5, 2) * 3 = (15, 6)
(0, 0, 0, 0)
向量的模
=(3,4) 的大小是多少？
根据勾股定理，的大小=
=5
我们用下面的符号(向量两边加双竖线)，来表示向量的长度，我们称该长度为向量的模
在三维空间中，
，我们首先知道
，则
，则最终得到
N维向量同理 
我们来扩展Attribute接口，增加一个求模的方法
public interface Attribute<T extends Number> { /** * 向量的长度 * @return */ int len(); /** * 获取向量的元素 * @param index 索引 * @return */ T get(int index); /** * 向量相加 * @param another * @return */ Attribute<T> add(Attribute<T> another,Addr<T> addr); /** * 向量乘以一个数 * @param k * @param addr * @return */ Attribute<T> mul(T k,Addr<T> addr); /** * 向量的模 * @return */ double norm(Addr<T> addr); }
Addr中增加一个平方和开方的方法
public interface Addr<T extends Number> { T zero(); T add(T lhs,T rhs); T mul(T k,T hs); T square(T hs); double prescription(T hs); }
向量实现类
public class Vector<T extends Number> implements Attribute<T> { private T[] values; public Vector(T[] values) { this.values = values; } @SuppressWarnings("unchecked") public static <T extends Number> Attribute<T> zero(int dim,Addr<T> addr) { T[] values = (T[])new Number[dim]; Arrays.fill(values,addr.zero()); return new Vector<>(values); } @Override public String toString() { return "Vector{" + "values=" + Arrays.toString(values) + '}'; } @Override public int len() { return values.length; } @Override public T get(int index) { if (index >= len() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } return values[index]; } @Override @SuppressWarnings("unchecked") public Attribute<T> add(Attribute<T> another,Addr<T> addr) { if (this.len() != another.len()) { throw new IllegalArgumentException("加法错误，两个向量的维度必须相等"); } T[] newValues = (T[])new Number[this.len()]; for (int i = 0; i < this.len(); i++) { List<T> list = new ArrayList<>(); list.add(this.get(i)); list.add(another.get(i)); newValues[i] = sumValue(list, addr); } this.values = newValues; return this; } @Override @SuppressWarnings("unchecked") public Attribute<T> mul(T k, Addr<T> addr) { T[] newValues = (T[])new Number[this.len()]; for (int i = 0; i < this.len(); i++) { newValues[i] = addr.mul(k,this.get(i)); } this.values = newValues; return this; } @Override public double norm(Addr<T> addr) { T total = addr.zero(); for (int i = 0; i < this.len(); i++) { total = addr.add(total,addr.square(this.get(i))); } return addr.prescription(total); } private T sumValue(List<T> list, Addr<T> addr) { T total = addr.zero(); for (T n : list) { total = addr.add(total,n); } return total; } public static void main(String[] args) { Attribute<Integer> vector = new Vector<>(new Integer[]{5,2}); System.out.println(vector.toString()); System.out.println(vector.len()); System.out.println(vector.get(0)); Attribute<Integer> vector1 = new Vector<>(new Integer[]{6,3}); Addr<Integer> addr = new Addr<Integer>() { @Override public Integer zero() { return 0; } @Override public Integer add(Integer lhs, Integer rhs) { return lhs + rhs; } @Override public Integer mul(Integer k, Integer hs) { return k * hs; } @Override public Integer square(Integer hs) { return hs * hs; } @Override public double prescription(Integer hs) { return Math.sqrt(hs); } }; vector.add(vector1, addr); System.out.println(vector); vector.mul(3,addr); System.out.println(vector); System.out.println(Vector.zero(4,addr)); System.out.println(vector.norm(addr)); } }
运行结果
Vector{values=[5, 2]}
2
5
Vector{values=[11, 5]}
Vector{values=[33, 15]}
Vector{values=[0, 0, 0, 0]}
36.24913792078372
Python代码
import math class Vector: def __init__(self, lst): self._values = lst @classmethod def zero(cls,dim): # 返回一个dim维的零向量 return cls([0] * dim) def norm(self): # 返回向量的模 return math.sqrt(sum(e**2 for e in self)) def __add__(self, another): # 向量加法，返回结果向量 assert len(self) == len(another),\ "加法错误，两个向量的维度必须相等" return Vector([a + b for a, b in zip(self, another)]) def __mul__(self, k): # 返回数量乘法的结果向量 return Vector([k * e for e in self]) def __iter__(self): # 返回向量的迭代器 return self._values.__iter__() def __getitem__(self, index): # 取向量的第index个元素 return self._values[index] def __len__(self): # 返回向量的长度 return len(self._values) def __repr__(self): return "Vector({})".format(self._values) def __str__(self): return "({})".format(", ".join(str(e) for e in self._values))
from playLA.Vector import Vector if __name__ == "__main__": vec = Vector([5, 2]) print(vec) print(len(vec)) print("vec[0] = {}, vec[1] = {}".format(vec[0], vec[1])) vec2 = Vector([6, 3]) print("{} + {} = {}".format(vec, vec2, vec + vec2)) print("{} * {} = {}".format(vec, 3, vec * 3)) zero2 = Vector.zero(4) print(zero2) print("norm({}) = {}".format(vec,vec.norm()))
运行结果
(5, 2)
2
vec[0] = 5, vec[1] = 2
(5, 2) + (6, 3) = (11, 5)
(5, 2) * 3 = (15, 6)
(0, 0, 0, 0)
norm((5, 2)) = 5.385164807134504
单位向量(unit vector)
就是向量的模为1的向量，所以每一个向量都有一个单位向量
对于一个N维向量来说，它的单位向量就是
，意思就是说对向量每一个维度的值去除以该向量的模，即为该向量的单位向量，这里我们需要记得，单位向量的模为1
，一般来说，单位向量只表示方向。我们把根据求出
的过程称为归一化，规范化(normalize)
图中蓝色的部分就是(3,4)这个向量的单位向量
单位向量有无数个，比如在二维空间中
我们以原点为圆心，1为半径画一个圆，那么这个圆中的所有点都是一个单位向量。
在二维空间中有两个特殊的单位向量,如下图中蓝线所示
=(1,0) 
=(0,1)，像这种只由0、1组成的单位向量，又叫标准单位向量(Standard Unit Vector)，标准单位向量指向坐标轴的正方向。
在三维空间中有三个标准单位向量
=(1,0,0) =(0,1,0) 
=(0,0,1)
N维空间中，有n个标准单位向量
=(1,0,...,0) =(0,1,...,0) ... 
=(0,0,...,1)
现在我们来扩展Attribute接口，新增一个单位向量的方法。
public interface Attribute<T extends Number> { /** * 向量的长度 * @return */ int len(); /** * 获取向量的元素 * @param index 索引 * @return */ T get(int index); /** * 向量相加 * @param another * @return */ Attribute<T> add(Attribute<T> another,Addr<T> addr); /** * 向量乘以一个数 * @param k * @param addr * @return */ Attribute<T> mul(T k,Addr<T> addr); /** * 向量的模 * @return */ double norm(Addr<T> addr); /** * 单位向量 * @return */ Attribute<Double> normalize(Addr<T> addr); }
addr中增加一个除法方法
public interface Addr<T extends Number> { T zero(); T add(T lhs,T rhs); T mul(T k,T hs); T square(T hs); double prescription(T hs); double div(T hs,double nhs); }
向量实现类
public class Vector<T extends Number> implements Attribute<T> { private T[] values; public Vector(T[] values) { this.values = values; } @SuppressWarnings("unchecked") public static <T extends Number> Attribute<T> zero(int dim,Addr<T> addr) { T[] values = (T[])new Number[dim]; Arrays.fill(values,addr.zero()); return new Vector<>(values); } @Override public String toString() { return "Vector{" + "values=" + Arrays.toString(values) + '}'; } @Override public int len() { return values.length; } @Override public T get(int index) { if (index >= len() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } return values[index]; } @Override @SuppressWarnings("unchecked") public Attribute<T> add(Attribute<T> another,Addr<T> addr) { if (this.len() != another.len()) { throw new IllegalArgumentException("加法错误，两个向量的维度必须相等"); } T[] newValues = (T[])new Number[this.len()]; for (int i = 0; i < this.len(); i++) { List<T> list = new ArrayList<>(); list.add(this.get(i)); list.add(another.get(i)); newValues[i] = sumValue(list, addr); } this.values = newValues; return this; } @Override @SuppressWarnings("unchecked") public Attribute<T> mul(T k, Addr<T> addr) { T[] newValues = (T[])new Number[this.len()]; for (int i = 0; i < this.len(); i++) { newValues[i] = addr.mul(k,this.get(i)); } this.values = newValues; return this; } @Override public double norm(Addr<T> addr) { T total = addr.zero(); for (int i = 0; i < this.len(); i++) { total = addr.add(total,addr.square(this.get(i))); } return addr.prescription(total); } @Override public Attribute<Double> normalize(Addr<T> addr) { Double[] newValues = new Double[this.len()]; double norm = this.norm(addr); for (int i = 0; i < this.len(); i++) { newValues[i] = addr.div(this.get(i),norm); } return new Vector<>(newValues); } private T sumValue(List<T> list, Addr<T> addr) { T total = addr.zero(); for (T n : list) { total = addr.add(total,n); } return total; } public static void main(String[] args) { Attribute<Integer> vector = new Vector<>(new Integer[]{5,2}); System.out.println(vector.toString()); System.out.println(vector.len()); System.out.println(vector.get(0)); Attribute<Integer> vector1 = new Vector<>(new Integer[]{6,3}); Addr<Integer> addr = new Addr<Integer>() { @Override public Integer zero() { return 0; } @Override public Integer add(Integer lhs, Integer rhs) { return lhs + rhs; } @Override public Integer mul(Integer k, Integer hs) { return k * hs; } @Override public Integer square(Integer hs) { return hs * hs; } @Override public double prescription(Integer hs) { return Math.sqrt(hs); } @Override public double div(Integer hs, double nhs) { return hs / nhs; } }; Addr<Double> addr1 = new Addr<Double>() { @Override public Double zero() { return 0.0; } @Override public Double add(Double lhs, Double rhs) { return lhs + rhs; } @Override public Double mul(Double k, Double hs) { return k * hs; } @Override public Double square(Double hs) { return hs * hs; } @Override public double prescription(Double hs) { return Math.sqrt(hs); } @Override public double div(Double hs, double nhs) { return hs / nhs; } }; vector.add(vector1, addr); System.out.println(vector); vector.mul(3,addr); System.out.println(vector); System.out.println(Vector.zero(4,addr)); System.out.println(vector.norm(addr)); Attribute<Double> unitVector = vector.normalize(addr); System.out.println(unitVector); System.out.println(unitVector.norm(addr1)); } }
运行结果
Vector{values=[5, 2]}
2
5
Vector{values=[11, 5]}
Vector{values=[33, 15]}
Vector{values=[0, 0, 0, 0]}
36.24913792078372
Vector{values=[0.9103664774626047, 0.413802944301184]}
0.9999999999999999
Python代码
import math class Vector: def __init__(self, lst): self._values = lst @classmethod def zero(cls,dim): # 返回一个dim维的零向量 return cls([0] * dim) def norm(self): # 返回向量的模 return math.sqrt(sum(e**2 for e in self)) def normalize(self): # 返回向量的单位向量 return Vector([e / self.norm() for e in self]) def __add__(self, another): # 向量加法，返回结果向量 assert len(self) == len(another),\ "加法错误，两个向量的维度必须相等" return Vector([a + b for a, b in zip(self, another)]) def __mul__(self, k): # 返回数量乘法的结果向量 return Vector([k * e for e in self]) def __iter__(self): # 返回向量的迭代器 return self._values.__iter__() def __getitem__(self, index): # 取向量的第index个元素 return self._values[index] def __len__(self): # 返回向量的长度 return len(self._values) def __repr__(self): return "Vector({})".format(self._values) def __str__(self): return "({})".format(", ".join(str(e) for e in self._values))
from playLA.Vector import Vector if __name__ == "__main__": vec = Vector([5, 2]) print(vec) print(len(vec)) print("vec[0] = {}, vec[1] = {}".format(vec[0], vec[1])) vec2 = Vector([6, 3]) print("{} + {} = {}".format(vec, vec2, vec + vec2)) print("{} * {} = {}".format(vec, 3, vec * 3)) zero2 = Vector.zero(4) print(zero2) print("norm({}) = {}".format(vec, vec.norm())) print("normalize {} is {}".format(vec, vec.normalize())) print(vec.normalize().norm())
运行结果
(5, 2)
2
vec[0] = 5, vec[1] = 2
(5, 2) + (6, 3) = (11, 5)
(5, 2) * 3 = (15, 6)
(0, 0, 0, 0)
norm((5, 2)) = 5.385164807134504
normalize (5, 2) is (0.9284766908852594, 0.3713906763541037)
1.0
向量的点乘
两个向量相乘
我们先给出它的计算公式，然来再来说明
由结果，我们可以看出，两个向量相乘，其结果是一个数(标量)。更严格的说法，这是两个向量的点乘，又或者叫两个向量的内积。
实际上，这个结果又
，即这两个向量的模乘以它们夹角的余弦。
我们在二维空间中来说明这个结论
在二维空间中，上述式子就转化成如下
我们假设它是这个样子的，由之前所知，是
与绿色向量的对角线，满足平行四边形法则，所以是+绿色向量，所以绿色向量=
，根据余弦定理，就有了
将式子调整一下就有了
根据模的计算方法，等式的右边展开，就有了
如此，我们就得到了上面的结论
在N维向量中，该等式依然成立
- 几何意义
我们来看一下的几何意义
我们以向方向上做一个投影，它的长度就是
。这样一来，我们就可以把等式的右边以这样来看待
，也就是说向量的点乘就是在方向上投影的大小乘以的大小，即让两个向量都指向同一个方向，在同一个方向上大小的乘积。
同理，我们也可以用在的方向上做投影，为
，换句话说，等式的右边以这样来看待
我们也可以将、分别投影到x轴和y轴上
排列组合后会得到4种情况
只有同方向上的投影相乘才是有意义的，而垂直方向上，它们彼此没有一点贡献，则为0.所以最终就只有
和
。
我们对Attribute接口进行一个扩展，新增一个点乘的方法。
public interface Attribute<T extends Number> { /** * 向量的长度 * @return */ int len(); /** * 获取向量的元素 * @param index 索引 * @return */ T get(int index); /** * 向量相加 * @param another * @return */ Attribute<T> add(Attribute<T> another,Addr<T> addr); /** * 向量乘以一个数 * @param k * @param addr * @return */ Attribute<T> mul(T k,Addr<T> addr); /** * 向量的模 * @return */ double norm(Addr<T> addr); /** * 单位向量 * @return */ Attribute<Double> normalize(Addr<T> addr); /** * 向量的点乘 * @param another * @param addr * @return */ T dot(Attribute<T> another,Addr<T> addr); }
向量实现类
public class Vector<T extends Number> implements Attribute<T> { private T[] values; public Vector(T[] values) { this.values = values; } @SuppressWarnings("unchecked") public static <T extends Number> Attribute<T> zero(int dim,Addr<T> addr) { T[] values = (T[])new Number[dim]; Arrays.fill(values,addr.zero()); return new Vector<>(values); } @Override public String toString() { return "Vector{" + "values=" + Arrays.toString(values) + '}'; } @Override public int len() { return values.length; } @Override public T get(int index) { if (index >= len() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } return values[index]; } @Override @SuppressWarnings("unchecked") public Attribute<T> add(Attribute<T> another,Addr<T> addr) { if (this.len() != another.len()) { throw new IllegalArgumentException("加法错误，两个向量的维度必须相等"); } T[] newValues = (T[])new Number[this.len()]; for (int i = 0; i < this.len(); i++) { List<T> list = new ArrayList<>(); list.add(this.get(i)); list.add(another.get(i)); newValues[i] = sumValue(list, addr); } this.values = newValues; return this; } @Override @SuppressWarnings("unchecked") public Attribute<T> mul(T k, Addr<T> addr) { T[] newValues = (T[])new Number[this.len()]; for (int i = 0; i < this.len(); i++) { newValues[i] = addr.mul(k,this.get(i)); } this.values = newValues; return this; } @Override public double norm(Addr<T> addr) { T total = addr.zero(); for (int i = 0; i < this.len(); i++) { total = addr.add(total,addr.square(this.get(i))); } return addr.prescription(total); } @Override public Attribute<Double> normalize(Addr<T> addr) { Double[] newValues = new Double[this.len()]; double norm = this.norm(addr); for (int i = 0; i < this.len(); i++) { newValues[i] = addr.div(this.get(i),norm); } return new Vector<>(newValues); } @Override public T dot(Attribute<T> another, Addr<T> addr) { if (this.len() != another.len()) { throw new IllegalArgumentException("点乘错误，两个向量的维度必须相等"); } T total = addr.zero(); for (int i = 0; i < this.len(); i++) { total = addr.add(total,addr.mul(this.get(i),another.get(i))); } return total; } private T sumValue(List<T> list, Addr<T> addr) { T total = addr.zero(); for (T n : list) { total = addr.add(total,n); } return total; } public static void main(String[] args) { Attribute<Integer> vector = new Vector<>(new Integer[]{5,2}); System.out.println(vector.toString()); System.out.println(vector.len()); System.out.println(vector.get(0)); Attribute<Integer> vector1 = new Vector<>(new Integer[]{6,3}); Addr<Integer> addr = new Addr<Integer>() { @Override public Integer zero() { return 0; } @Override public Integer add(Integer lhs, Integer rhs) { return lhs + rhs; } @Override public Integer mul(Integer k, Integer hs) { return k * hs; } @Override public Integer square(Integer hs) { return hs * hs; } @Override public double prescription(Integer hs) { return Math.sqrt(hs); } @Override public double div(Integer hs, double nhs) { return hs / nhs; } }; Addr<Double> addr1 = new Addr<Double>() { @Override public Double zero() { return 0.0; } @Override public Double add(Double lhs, Double rhs) { return lhs + rhs; } @Override public Double mul(Double k, Double hs) { return k * hs; } @Override public Double square(Double hs) { return hs * hs; } @Override public double prescription(Double hs) { return Math.sqrt(hs); } @Override public double div(Double hs, double nhs) { return hs / nhs; } }; //(5,2)*(6,3) System.out.println(vector.dot(vector1,addr)); //(5,2)+(6,3) vector.add(vector1, addr); System.out.println(vector); vector.mul(3,addr); System.out.println(vector); System.out.println(Vector.zero(4,addr)); System.out.println(vector.norm(addr)); Attribute<Double> unitVector = vector.normalize(addr); System.out.println(unitVector); System.out.println(unitVector.norm(addr1)); } }
运行结果
Vector{values=[5, 2]}
2
5
36
Vector{values=[11, 5]}
Vector{values=[33, 15]}
Vector{values=[0, 0, 0, 0]}
36.24913792078372
Vector{values=[0.9103664774626047, 0.413802944301184]}
0.9999999999999999
Python代码
import math class Vector: def __init__(self, lst): self._values = lst @classmethod def zero(cls,dim): # 返回一个dim维的零向量 return cls([0] * dim) def norm(self): # 返回向量的模 return math.sqrt(sum(e**2 for e in self)) def normalize(self): # 返回向量的单位向量 return Vector([e / self.norm() for e in self]) def __add__(self, another): # 向量加法，返回结果向量 assert len(self) == len(another),\ "加法错误，两个向量的维度必须相等" return Vector([a + b for a, b in zip(self, another)]) def dot(self, another): # 向量点乘，返回结果标量 assert len(self) == len(another), \ "点乘错误，两个向量的维度必须相等" return sum(a * b for a, b in zip(self, another)) def __mul__(self, k): # 返回数量乘法的结果向量 return Vector([k * e for e in self]) def __iter__(self): # 返回向量的迭代器 return self._values.__iter__() def __getitem__(self, index): # 取向量的第index个元素 return self._values[index] def __len__(self): # 返回向量的长度 return len(self._values) def __repr__(self): return "Vector({})".format(self._values) def __str__(self): return "({})".format(", ".join(str(e) for e in self._values))
from playLA.Vector import Vector if __name__ == "__main__": vec = Vector([5, 2]) print(vec) print(len(vec)) print("vec[0] = {}, vec[1] = {}".format(vec[0], vec[1])) vec2 = Vector([6, 3]) print("{} + {} = {}".format(vec, vec2, vec + vec2)) print("{} * {} = {}".format(vec, 3, vec * 3)) zero2 = Vector.zero(4) print(zero2) print("norm({}) = {}".format(vec, vec.norm())) print("normalize {} is {}".format(vec, vec.normalize())) print(vec.normalize().norm()) print(vec.dot(vec2))
运行结果
(5, 2)
2
vec[0] = 5, vec[1] = 2
(5, 2) + (6, 3) = (11, 5)
(5, 2) * 3 = (15, 6)
(0, 0, 0, 0)
norm((5, 2)) = 5.385164807134504
normalize (5, 2) is (0.9284766908852594, 0.3713906763541037)
1.0
36
向量点乘的应用
由该公式，我们可以得出两个向量之间的夹角的余弦等于两个向量的点乘除以它们模的积
如果
，
。
反之，则两个向量是垂直的。知道两个向量是否是垂直的，是一件非常重要的事情。
如果
，两个向量夹角为锐角。
如果
，两个向量夹角为钝角。
计算两个向量的夹角
Attribute新增两个向量的夹角的方法
public interface Attribute<T extends Number> { /** * 向量的长度 * @return */ int len(); /** * 获取向量的元素 * @param index 索引 * @return */ T get(int index); /** * 向量相加 * @param another * @return */ Attribute<T> add(Attribute<T> another,Addr<T> addr); /** * 向量乘以一个数 * @param k * @param addr * @return */ Attribute<T> mul(T k,Addr<T> addr); /** * 向量的模 * @return */ double norm(Addr<T> addr); /** * 单位向量 * @return */ Attribute<Double> normalize(Addr<T> addr); /** * 向量的点乘 * @param another * @param addr * @return */ T dot(Attribute<T> another,Addr<T> addr); /** * 向量的夹角 * @param another * @param addr * @return */ double degree(Attribute<T> another,Addr<T> addr); }
Addr新增反余弦和度数解析方法
public interface Addr<T extends Number> { T zero(); T add(T lhs,T rhs); T mul(T k,T hs); T square(T hs); double prescription(T hs); double div(T hs,double nhs); double acos(double hs); double toDegree(double hs); }
向量实现类
public class Vector<T extends Number> implements Attribute<T> { private T[] values; public Vector(T[] values) { this.values = values; } @SuppressWarnings("unchecked") public static <T extends Number> Attribute<T> zero(int dim,Addr<T> addr) { T[] values = (T[])new Number[dim]; Arrays.fill(values,addr.zero()); return new Vector<>(values); } @Override public String toString() { return "Vector{" + "values=" + Arrays.toString(values) + '}'; } @Override public int len() { return values.length; } @Override public T get(int index) { if (index >= len() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } return values[index]; } @Override @SuppressWarnings("unchecked") public Attribute<T> add(Attribute<T> another,Addr<T> addr) { if (this.len() != another.len()) { throw new IllegalArgumentException("加法错误，两个向量的维度必须相等"); } T[] newValues = (T[])new Number[this.len()]; for (int i = 0; i < this.len(); i++) { List<T> list = new ArrayList<>(); list.add(this.get(i)); list.add(another.get(i)); newValues[i] = sumValue(list, addr); } this.values = newValues; return this; } @Override @SuppressWarnings("unchecked") public Attribute<T> mul(T k, Addr<T> addr) { T[] newValues = (T[])new Number[this.len()]; for (int i = 0; i < this.len(); i++) { newValues[i] = addr.mul(k,this.get(i)); } this.values = newValues; return this; } @Override public double norm(Addr<T> addr) { T total = addr.zero(); for (int i = 0; i < this.len(); i++) { total = addr.add(total,addr.square(this.get(i))); } return addr.prescription(total); } @Override public Attribute<Double> normalize(Addr<T> addr) { Double[] newValues = new Double[this.len()]; double norm = this.norm(addr); for (int i = 0; i < this.len(); i++) { newValues[i] = addr.div(this.get(i),norm); } return new Vector<>(newValues); } @Override public T dot(Attribute<T> another, Addr<T> addr) { if (this.len() != another.len()) { throw new IllegalArgumentException("点乘错误，两个向量的维度必须相等"); } T total = addr.zero(); for (int i = 0; i < this.len(); i++) { total = addr.add(total,addr.mul(this.get(i),another.get(i))); } return total; } @Override public double degree(Attribute<T> another, Addr<T> addr) { T dot = this.dot(another,addr); double selfNorm = this.norm(addr); double anotherNorm = another.norm(addr); double productNorm = selfNorm * anotherNorm; return addr.toDegree(addr.acos(addr.div(dot,productNorm))); } private T sumValue(List<T> list, Addr<T> addr) { T total = addr.zero(); for (T n : list) { total = addr.add(total,n); } return total; } public static void main(String[] args) { Attribute<Integer> vector = new Vector<>(new Integer[]{5,2}); System.out.println(vector.toString()); System.out.println(vector.len()); System.out.println(vector.get(0)); Attribute<Integer> vector1 = new Vector<>(new Integer[]{6,3}); Addr<Integer> addr = new Addr<Integer>() { @Override public Integer zero() { return 0; } @Override public Integer add(Integer lhs, Integer rhs) { return lhs + rhs; } @Override public Integer mul(Integer k, Integer hs) { return k * hs; } @Override public Integer square(Integer hs) { return hs * hs; } @Override public double prescription(Integer hs) { return Math.sqrt(hs); } @Override public double div(Integer hs, double nhs) { return hs / nhs; } @Override public double acos(double hs) { return Math.acos(hs); } @Override public double toDegree(double hs) { return Math.toDegrees(hs); } }; Addr<Double> addr1 = new Addr<Double>() { @Override public Double zero() { return 0.0; } @Override public Double add(Double lhs, Double rhs) { return lhs + rhs; } @Override public Double mul(Double k, Double hs) { return k * hs; } @Override public Double square(Double hs) { return hs * hs; } @Override public double prescription(Double hs) { return Math.sqrt(hs); } @Override public double div(Double hs, double nhs) { return hs / nhs; } @Override public double acos(double hs) { return Math.acos(hs); } @Override public double toDegree(double hs) { return Math.toDegrees(hs); } }; //(5,2)*(6,3) System.out.println(vector.dot(vector1,addr)); //(5,2),(6,3)的夹角 System.out.println(vector.degree(vector1,addr)); //(5,2)+(6,3) vector.add(vector1, addr); System.out.println(vector); vector.mul(3,addr); System.out.println(vector); System.out.println(Vector.zero(4,addr)); System.out.println(vector.norm(addr)); Attribute<Double> unitVector = vector.normalize(addr); System.out.println(unitVector); System.out.println(unitVector.norm(addr1)); } }
运行结果
Vector{values=[5, 2]}
2
5
36
4.763641690726143
Vector{values=[11, 5]}
Vector{values=[33, 15]}
Vector{values=[0, 0, 0, 0]}
36.24913792078372
Vector{values=[0.9103664774626047, 0.413802944301184]}
0.9999999999999999
Python代码
import math class Vector: def __init__(self, lst): self._values = lst @classmethod def zero(cls,dim): # 返回一个dim维的零向量 return cls([0] * dim) def degree(self, another): # 返回向量的夹角 return math.degrees(math.acos(self.dot(another) / (self.norm() * another.norm()))) def norm(self): # 返回向量的模 return math.sqrt(sum(e**2 for e in self)) def normalize(self): # 返回向量的单位向量 return Vector([e / self.norm() for e in self]) def __add__(self, another): # 向量加法，返回结果向量 assert len(self) == len(another),\ "加法错误，两个向量的维度必须相等" return Vector([a + b for a, b in zip(self, another)]) def dot(self, another): # 向量点乘，返回结果标量 assert len(self) == len(another), \ "点乘错误，两个向量的维度必须相等" return sum(a * b for a, b in zip(self, another)) def __mul__(self, k): # 返回数量乘法的结果向量 return Vector([k * e for e in self]) def __iter__(self): # 返回向量的迭代器 return self._values.__iter__() def __getitem__(self, index): # 取向量的第index个元素 return self._values[index] def __len__(self): # 返回向量的长度 return len(self._values) def __repr__(self): return "Vector({})".format(self._values) def __str__(self): return "({})".format(", ".join(str(e) for e in self._values))
from playLA.Vector import Vector if __name__ == "__main__": vec = Vector([5, 2]) print(vec) print(len(vec)) print("vec[0] = {}, vec[1] = {}".format(vec[0], vec[1])) vec2 = Vector([6, 3]) print("{} + {} = {}".format(vec, vec2, vec + vec2)) print("{} * {} = {}".format(vec, 3, vec * 3)) zero2 = Vector.zero(4) print(zero2) print("norm({}) = {}".format(vec, vec.norm())) print("normalize {} is {}".format(vec, vec.normalize())) print(vec.normalize().norm()) print(vec.dot(vec2)) print(vec.degree(vec2))
运行结果
(5, 2)
2
vec[0] = 5, vec[1] = 2
(5, 2) + (6, 3) = (11, 5)
(5, 2) * 3 = (15, 6)
(0, 0, 0, 0)
norm((5, 2)) = 5.385164807134504
normalize (5, 2) is (0.9284766908852594, 0.3713906763541037)
1.0
36
4.763641690726144
在N维向量的标准单位向量一定是两两互相垂直的
在二维空间中
=(1,0) =(0,1) 
=0
三维空间中
=(1,0,0) =(0,1,0) =(0,0,1)
=0 
=0 
=0
判断两个向量的相似程度(推荐系统)
在推荐系统中针对一组物品进行推荐，这组物品有多个维度的属性，我们把这组物品的属性组合成一个多维空间的向量。如果你在一个系统中留下了足迹，那么系统就会根据你的喜好，把相似的属性的物品(多维空间向量)推荐给你。当然这个是推荐系统最基本的理论——余弦相似度，现在成熟的推荐系统要比这个复杂的多，但是这个是最基本的理论基础。
两个物品的多维属性形成的向量的夹角如果是锐角，我们就说这两个物品相似
相似
两个物品的多维属性形成的向量的夹角如果是直角，我们就说这两个物品无关
无关
两个物品的多维属性形成的向量的夹角如果是钝角，我们就说这两个物品背离
背离
在几何计算中，我们想知道一个向量在另一个向量中的投影的坐标
我们知道投影到上的长度为
，定义为d
则
之前我们知道单位向量只表示方向，它的模为1，则的单位向量为
，则投影点的坐标为
扩展Attribute接口，新增一个向量映射的方法
public interface Attribute<T extends Number> { /** * 向量的长度 * @return */ int len(); /** * 获取向量的元素 * @param index 索引 * @return */ T get(int index); /** * 向量相加 * @param another * @return */ Attribute<T> add(Attribute<T> another,Addr<T> addr); /** * 向量乘以一个数 * @param k * @param addr * @return */ Attribute<T> mul(T k,Addr<T> addr); /** * 向量的模 * @return */ double norm(Addr<T> addr); /** * 单位向量 * @return */ Attribute<Double> normalize(Addr<T> addr); /** * 向量的点乘 * @param another * @param addr * @return */ T dot(Attribute<T> another,Addr<T> addr); /** * 向量的夹角 * @param another * @param addr * @return */ double degree(Attribute<T> another,Addr<T> addr); /** * 其他向量在当前向量的投影 * @param another * @param addr * @return */ Attribute<Double> project(Attribute<T> another,Addr<T> addr); }
向量实现类
public class Vector<T extends Number> implements Attribute<T> { private T[] values; public Vector(T[] values) { this.values = values; } @SuppressWarnings("unchecked") public static <T extends Number> Attribute<T> zero(int dim,Addr<T> addr) { T[] values = (T[])new Number[dim]; Arrays.fill(values,addr.zero()); return new Vector<>(values); } @Override public String toString() { return "Vector{" + "values=" + Arrays.toString(values) + '}'; } @Override public int len() { return values.length; } @Override public T get(int index) { if (index >= len() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } return values[index]; } @Override @SuppressWarnings("unchecked") public Attribute<T> add(Attribute<T> another,Addr<T> addr) { if (this.len() != another.len()) { throw new IllegalArgumentException("加法错误，两个向量的维度必须相等"); } T[] newValues = (T[])new Number[this.len()]; for (int i = 0; i < this.len(); i++) { List<T> list = new ArrayList<>(); list.add(this.get(i)); list.add(another.get(i)); newValues[i] = sumValue(list, addr); } this.values = newValues; return this; } @Override @SuppressWarnings("unchecked") public Attribute<T> mul(T k, Addr<T> addr) { T[] newValues = (T[])new Number[this.len()]; for (int i = 0; i < this.len(); i++) { newValues[i] = addr.mul(k,this.get(i)); } this.values = newValues; return this; } @Override public double norm(Addr<T> addr) { T total = addr.zero(); for (int i = 0; i < this.len(); i++) { total = addr.add(total,addr.square(this.get(i))); } return addr.prescription(total); } @Override public Attribute<Double> normalize(Addr<T> addr) { Double[] newValues = new Double[this.len()]; double norm = this.norm(addr); for (int i = 0; i < this.len(); i++) { newValues[i] = addr.div(this.get(i),norm); } return new Vector<>(newValues); } @Override public T dot(Attribute<T> another, Addr<T> addr) { if (this.len() != another.len()) { throw new IllegalArgumentException("点乘错误，两个向量的维度必须相等"); } T total = addr.zero(); for (int i = 0; i < this.len(); i++) { total = addr.add(total,addr.mul(this.get(i),another.get(i))); } return total; } @Override public double degree(Attribute<T> another, Addr<T> addr) { T dot = this.dot(another,addr); double selfNorm = this.norm(addr); double anotherNorm = another.norm(addr); double productNorm = selfNorm * anotherNorm; return addr.toDegree(addr.acos(addr.div(dot,productNorm))); } @Override public Attribute<Double> project(Attribute<T> another, Addr<T> addr) { T dot = this.dot(another,addr); double selfNorm = this.norm(addr); Attribute<Double> unitVector = this.normalize(addr); double d = addr.div(dot, selfNorm); return mul(d,unitVector); } private Attribute<Double> mul(double d,Attribute<Double> v) { Double[] values = new Double[v.len()]; for (int i = 0; i < v.len(); i++) { values[i] = d * v.get(i); } return new Vector<>(values); } private T sumValue(List<T> list, Addr<T> addr) { T total = addr.zero(); for (T n : list) { total = addr.add(total,n); } return total; } public static void main(String[] args) { Attribute<Integer> vector = new Vector<>(new Integer[]{5,2}); System.out.println(vector.toString()); System.out.println(vector.len()); System.out.println(vector.get(0)); Attribute<Integer> vector1 = new Vector<>(new Integer[]{6,3}); Addr<Integer> addr = new Addr<Integer>() { @Override public Integer zero() { return 0; } @Override public Integer add(Integer lhs, Integer rhs) { return lhs + rhs; } @Override public Integer mul(Integer k, Integer hs) { return k * hs; } @Override public Integer square(Integer hs) { return hs * hs; } @Override public double prescription(Integer hs) { return Math.sqrt(hs); } @Override public double div(Integer hs, double nhs) { return hs / nhs; } @Override public double acos(double hs) { return Math.acos(hs); } @Override public double toDegree(double hs) { return Math.toDegrees(hs); } }; Addr<Double> addr1 = new Addr<Double>() { @Override public Double zero() { return 0.0; } @Override public Double add(Double lhs, Double rhs) { return lhs + rhs; } @Override public Double mul(Double k, Double hs) { return k * hs; } @Override public Double square(Double hs) { return hs * hs; } @Override public double prescription(Double hs) { return Math.sqrt(hs); } @Override public double div(Double hs, double nhs) { return hs / nhs; } @Override public double acos(double hs) { return Math.acos(hs); } @Override public double toDegree(double hs) { return Math.toDegrees(hs); } }; //(5,2)*(6,3) System.out.println(vector.dot(vector1,addr)); //(5,2),(6,3)的夹角 System.out.println(vector.degree(vector1,addr)); //(6,3)在(5,2)上的映射 System.out.println(vector.project(vector1,addr)); //(5,2)+(6,3) vector.add(vector1, addr); System.out.println(vector); vector.mul(3,addr); System.out.println(vector); System.out.println(Vector.zero(4,addr)); System.out.println(vector.norm(addr)); Attribute<Double> unitVector = vector.normalize(addr); System.out.println(unitVector); System.out.println(unitVector.norm(addr1)); } }
运行结果
Vector{values=[5, 2]}
2
5
36
4.763641690726143
Vector{values=[6.20689655172414, 2.4827586206896552]}
Vector{values=[11, 5]}
Vector{values=[33, 15]}
Vector{values=[0, 0, 0, 0]}
36.24913792078372
Vector{values=[0.9103664774626047, 0.413802944301184]}
0.9999999999999999
Python代码
import math class Vector: def __init__(self, lst): self._values = lst @classmethod def zero(cls,dim): # 返回一个dim维的零向量 return cls([0] * dim) def project(self, another): # 其他向量在当前向量的投影 return self.normalize() * (self.dot(another) / self.norm()) def degree(self, another): # 返回向量的夹角 return math.degrees(math.acos(self.dot(another) / (self.norm() * another.norm()))) def norm(self): # 返回向量的模 return math.sqrt(sum(e**2 for e in self)) def normalize(self): # 返回向量的单位向量 return Vector([e / self.norm() for e in self]) def __add__(self, another): # 向量加法，返回结果向量 assert len(self) == len(another),\ "加法错误，两个向量的维度必须相等" return Vector([a + b for a, b in zip(self, another)]) def dot(self, another): # 向量点乘，返回结果标量 assert len(self) == len(another), \ "点乘错误，两个向量的维度必须相等" return sum(a * b for a, b in zip(self, another)) def __mul__(self, k): # 返回数量乘法的结果向量 return Vector([k * e for e in self]) def __iter__(self): # 返回向量的迭代器 return self._values.__iter__() def __getitem__(self, index): # 取向量的第index个元素 return self._values[index] def __len__(self): # 返回向量的长度 return len(self._values) def __repr__(self): return "Vector({})".format(self._values) def __str__(self): return "({})".format(", ".join(str(e) for e in self._values))
from playLA.Vector import Vector if __name__ == "__main__": vec = Vector([5, 2]) print(vec) print(len(vec)) print("vec[0] = {}, vec[1] = {}".format(vec[0], vec[1])) vec2 = Vector([6, 3]) print("{} + {} = {}".format(vec, vec2, vec + vec2)) print("{} * {} = {}".format(vec, 3, vec * 3)) zero2 = Vector.zero(4) print(zero2) print("norm({}) = {}".format(vec, vec.norm())) print("normalize {} is {}".format(vec, vec.normalize())) print(vec.normalize().norm()) print(vec.dot(vec2)) print(vec.degree(vec2)) print(vec.project(vec2))
运行结果
(5, 2)
2
vec[0] = 5, vec[1] = 2
(5, 2) + (6, 3) = (11, 5)
(5, 2) * 3 = (15, 6)
(0, 0, 0, 0)
norm((5, 2)) = 5.385164807134504
normalize (5, 2) is (0.9284766908852594, 0.3713906763541037)
1.0
36
4.763641690726144
(6.20689655172414, 2.4827586206896552)
Numpy中向量的基本使用
import numpy as np if __name__ == "__main__": print(np.__version__) # 在Python中，列表可以任意更改元素类型 lst = [1, 2, 3] lst[0] = "早上好" print(lst) # 在numpy中则不可以 vec = np.array([1, 2, 3]) print(vec) # numpy中的向量可以修改元素内容 # vec[0] = 666 print(vec) # 创建零向量 print(np.zeros(4)) print(np.ones(4)) print(np.full(4, 666)) # 获取向量的长度 print("size =", vec.size) print("size =", len(vec)) # 获取向量中的元素 print(vec[0]) # 获取向量最后一个元素 print(vec[-1]) # 获取向量中的一段元素，返回一个新向量 print(vec[0: 2]) print(type(vec[0: 2])) vec2 = np.array([4, 5, 6]) # 向量相加 print("{} + {} = {}".format(vec, vec2, vec + vec2)) # 向量相减 print("{} - {} = {}".format(vec, vec2, vec - vec2)) # 向量的数量乘 print("{} * {} = {}".format(2, vec, 2 * vec)) print("{} * {} = {}".format(vec, 2, vec * 2)) # 向量与向量的*号乘法，本身没有数学意义 print("{} * {} = {}".format(vec, vec2, vec * vec2)) # 向量与向量的点乘 print("{}.dot({}) = {}".format(vec, vec2, vec.dot(vec2))) # 向量的模 print(np.linalg.norm(vec)) # numpy中没有实现单位向量，需要自己实现 print(vec / np.linalg.norm(vec)) # 单位向量的模 print(np.linalg.norm(vec / np.linalg.norm(vec)))
运行结果
1.18.1
['早上好', 2, 3]
[1 2 3]
[1 2 3]
[0. 0. 0. 0.]
[1. 1. 1. 1.]
[666 666 666 666]
size = 3
size = 3
1
3
[1 2]
<class 'numpy.ndarray'>
[1 2 3] + [4 5 6] = [5 7 9]
[1 2 3] - [4 5 6] = [-3 -3 -3]
2 * [1 2 3] = [2 4 6]
[1 2 3] * 2 = [2 4 6]
[1 2 3] * [4 5 6] = [ 4 10 18]
[1 2 3].dot([4 5 6]) = 32
3.7416573867739413
[0.26726124 0.53452248 0.80178373]
1.0
矩阵
像这样排成行列的一组数，我们称之为矩阵。向量对数的拓展，一个向量表示一组数。矩阵是对向量的拓展，一个矩阵表示一组向量。
根据看待的视角不同，我们可以把矩阵以行、列来划分成不同的向量
对于这样的一个矩阵，我们可以把它称为一个4 * 4矩阵，因为它行数为4，列数为4
对于这样一个矩阵，我们则称它为3 * 4矩阵，因为它行数为3，列数为4.
对于行数 = 列数的矩阵，有一个特殊的名字——方阵，方阵有很多特殊的性质，有很多特殊的矩阵是方阵。
我们通常用A来代表一个矩阵，而矩阵中的元素用
来表示，其中i表示第i行，j表示第j列。
矩阵在不同的应用中，可以表示不同的意义。比如这张成绩表就可以看成一个矩阵，它的行向量是一个一个学生的各个科目的成绩。而它的列向量则是某一个具体科目的所有学生成绩。
现在我们来实现矩阵类，先定义一个接口
public interface Ranks<T extends Number> { /** * 矩阵的行列数 * @return */ int[] shape(); /** * 矩阵的行数 * @return */ int rowNum(); /** * 矩阵的列数 * @return */ int colNum(); /** * 矩阵的元素个数 * @return */ int size(); /** * 获取矩阵中的一个元素 * @param pos * @return */ T get(int pos); /** * 获取矩阵中的一个元素 * @param pos * @return */ T get(int[] pos); /** * 获取一个行向量 * @param index * @return */ Attribute<T> rowVector(int index); /** * 获取一个列向量 * @param index * @return */ Attribute<T> colVector(int index); }
矩阵实现类
public class Matrix<T extends Number> implements Ranks<T>{ private T[][] values; public Matrix(T[][] values) { this.values = values; } @Override public String toString() { StringBuilder builder = new StringBuilder(); for (int i = 0; i < values.length - 1; i++) { builder.append(Arrays.toString(values[i]) + ","); } builder.append(Arrays.toString(values[values.length - 1])); return "Matrix{" + "values=" + "[" + builder.toString() + "]" + '}'; } @Override public int[] shape() { int[] res = new int[2]; res[0] = values.length; res[1] = values[0].length; return res; } @Override public int rowNum() { return shape()[0]; } @Override public int colNum() { return shape()[1]; } @Override public int size() { int row = rowNum(); int col = colNum(); return row * col; } @Override public T get(int pos) { if (pos >= size() || pos < 0) { throw new IllegalArgumentException("位置超出索引范围"); } int row = pos / colNum(); int col = pos % colNum(); return values[row][col]; } @Override public T get(int[] pos) { if (pos.length != 2) { throw new IllegalArgumentException("位置长度超出范围"); } if (pos[0] >= rowNum() || pos[0] < 0 || pos[1] >= colNum() || pos[1] < 0) { throw new IllegalArgumentException("位置超出索引范围"); } return values[pos[0]][pos[1]]; } @Override public Attribute<T> rowVector(int index) { if (index >= rowNum() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } return new Vector<>(values[index]); } @Override @SuppressWarnings("unchecked") public Attribute<T> colVector(int index) { if (index >= colNum() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } T[] values = (T[])new Number[rowNum()]; for (int i = 0; i < rowNum(); i++) { values[i] = this.values[i][index]; } return new Vector<>(values); } public static void main(String[] args) { Ranks<Integer> matrix = new Matrix<>(new Integer[][]{{1,2,3},{3,4,5}}); System.out.println(matrix); System.out.println(Arrays.toString(matrix.shape())); System.out.println(matrix.rowNum()); System.out.println(matrix.colNum()); System.out.println(matrix.size()); System.out.println(matrix.get(2)); System.out.println(matrix.get(new int[]{0,0})); System.out.println(matrix.rowVector(1)); System.out.println(matrix.colVector(0)); } }
运行结果
Matrix{values=[[1, 2, 3],[3, 4, 5]]}
[2, 3]
2
3
6
3
1
Vector{values=[3, 4, 5]}
Vector{values=[1, 3]}
Python代码
from .Vector import Vector class Matrix: def __init__(self, list2d): self._values = [row[:] for row in list2d] def row_vector(self, index): # 返回矩阵的第index个行向量 return Vector(self._values[index]) def col_vector(self, index): # 返回矩阵的第index个列向量 return Vector([row[index] for row in self._values]) def __getitem__(self, pos): # 返回矩阵pos位置的元素 r, c = pos return self._values[r][c] def size(self): # 返回矩阵中元素的个数 r, c = self.shape() return r * c def row_num(self): # 返回矩阵的行数 return self.shape()[0] __len__ = row_num def col_num(self): # 返回矩阵的列数 return self.shape()[1] def shape(self): # 返回矩阵的形状：(行数，列数) return len(self._values), len(self._values[0]) def __repr__(self): return "Matrix({})".format(self._values) __str__ = __repr__
from playLA.Matrix import Matrix if __name__ == "__main__": matrix = Matrix([[1, 2, 3], [3, 4, 5]]) print(matrix) print("matrix.shape = {}".format(matrix.shape())) print("matrix.size = {}".format(matrix.size())) print("len(matrix) = {}".format(len(matrix))) print("matrix[0][0] = {}".format(matrix[0, 0])) print(matrix.row_vector(1)) print(matrix.col_vector(0))
运行结果
Matrix([[1, 2, 3], [3, 4, 5]])
matrix.shape = (2, 3)
matrix.size = 6
len(matrix) = 2
matrix[0][0] = 1
(3, 4, 5)
(1, 3)
矩阵的基本运算
类比于向量的基本运算，矩阵的基本运算也包括矩阵的加法A + B，矩阵的数量乘法
对于矩阵A+B
则有
而矩阵的数量乘法
假设有一个矩阵如下
则它乘以2之后为
它们的几何意义为
P即为二维平面的一个三角形，而2P则是将该三角形进行了一次缩放。矩阵的运算在计算机图形上有很多的应用。
矩阵的基本运算性质
加法交换律 A + B = B + A
加法结合律 (A + B) + C = A + (B + C)
存在矩阵O,满足 A + O = A
存在矩阵-A，满足 A + (-A) = O,-A唯一，-A = - 1 * A
乘法结合律 (ck)A = c(kA)
乘法分配律 k(A + B) = kA + kB
(c + k)A = cA + kA
扩充Ranks接口
public interface Ranks<T extends Number> { /** * 矩阵的行列数 * @return */ int[] shape(); /** * 矩阵的行数 * @return */ int rowNum(); /** * 矩阵的列数 * @return */ int colNum(); /** * 矩阵的元素个数 * @return */ int size(); /** * 获取矩阵中的一个元素 * @param pos * @return */ T get(int pos); /** * 获取矩阵中的一个元素 * @param pos * @return */ T get(int[] pos); /** * 获取一个行向量 * @param index * @return */ Attribute<T> rowVector(int index); /** * 获取一个列向量 * @param index * @return */ Attribute<T> colVector(int index); /** * 矩阵相加 * @param another * @return */ Ranks<T> add(Ranks<T> another,Addr<T> addr); /** * 矩阵相减 * @param another * @param addr * @return */ Ranks<T> sub(Ranks<T> another,Addr<T> addr); /** * 矩阵乘以一个数 * @param k * @param addr * @return */ Ranks<T> mul(T k,Addr<T> addr); /** * 矩阵除以一个数 * @param k * @param addr * @return */ Ranks<Double> div(double k,Addr<T> addr); /** * 矩阵取正的结果 * @return */ Ranks<T> pos(Addr<T> addr,T one); /** * 矩阵取负的结果 * @return */ Ranks<T> neg(Addr<T> addr,T negOne); }
addr增加减法运算
public interface Addr<T extends Number> { T zero(); T add(T lhs,T rhs); T sub(T lhs,T rhs); T mul(T k,T hs); T square(T hs); double prescription(T hs); double div(T hs,double nhs); double acos(double hs); double toDegree(double hs); }
矩阵实现类
public class Matrix<T extends Number> implements Ranks<T>,Cloneable { private T[][] values; public Matrix(T[][] values) { this.values = values; } @Override public String toString() { StringBuilder builder = new StringBuilder(); for (int i = 0; i < values.length - 1; i++) { builder.append(Arrays.toString(values[i]) + ","); } builder.append(Arrays.toString(values[values.length - 1])); return "Matrix{" + "values=" + "[" + builder.toString() + "]" + '}'; } /** * 创建一个零矩阵 * @param row * @param col * @param addr * @param <T> * @return */ @SuppressWarnings("unchecked") public static <T extends Number> Ranks<T> zero(int row,int col,Addr<T> addr) { T[][] values = (T[][])new Number[row][col]; for (int i = 0; i < row; i++) { Arrays.fill(values[i],addr.zero()); } return new Matrix<>(values); } @Override public int[] shape() { int[] res = new int[2]; res[0] = values.length; res[1] = values[0].length; return res; } @Override public int rowNum() { return shape()[0]; } @Override public int colNum() { return shape()[1]; } @Override public int size() { int row = rowNum(); int col = colNum(); return row * col; } @Override public T get(int pos) { if (pos >= size() || pos < 0) { throw new IllegalArgumentException("位置超出索引范围"); } int row = pos / colNum(); int col = pos % colNum(); return values[row][col]; } @Override public T get(int[] pos) { if (pos.length != 2) { throw new IllegalArgumentException("位置长度超出范围"); } if (pos[0] >= rowNum() || pos[0] < 0 || pos[1] >= colNum() || pos[1] < 0) { throw new IllegalArgumentException("位置超出索引范围"); } return values[pos[0]][pos[1]]; } @Override public Attribute<T> rowVector(int index) { if (index >= rowNum() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } return new Vector<>(values[index]); } @Override @SuppressWarnings("unchecked") public Attribute<T> colVector(int index) { if (index >= colNum() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } T[] values = (T[])new Number[rowNum()]; for (int i = 0; i < rowNum(); i++) { values[i] = this.values[i][index]; } return new Vector<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> add(Ranks<T> another,Addr<T> addr) { if (!Arrays.equals(this.shape(),another.shape())) { throw new IllegalArgumentException("加法错误，两个矩阵的形状必须相同"); } T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.add(this.get(new int[]{i,j}),another.get(new int[]{i,j})); } } this.values = values; return this; } @Override @SuppressWarnings("unchecked") public Ranks<T> sub(Ranks<T> another, Addr<T> addr) { if (!Arrays.equals(this.shape(),another.shape())) { throw new IllegalArgumentException("减法错误，两个矩阵的形状必须相同"); } T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.sub(this.get(new int[]{i,j}),another.get(new int[]{i,j})); } } this.values = values; return this; } @Override @SuppressWarnings("unchecked") public Ranks<T> mul(T k, Addr<T> addr) { T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.mul(k,this.get(new int[]{i,j})); } } this.values = values; return this; } @Override public Ranks<Double> div(double k, Addr<T> addr) { Double[][] values = new Double[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.div(this.get(new int[]{i,j}),k); } } return new Matrix<>(values); } @Override public Ranks<T> pos(Addr<T> addr,T one) { return mul(one,addr); } @Override public Ranks<T> neg(Addr<T> addr,T negOne) { return mul(negOne,addr); } @Override public Matrix clone() throws CloneNotSupportedException { return (Matrix)super.clone(); } @SuppressWarnings("unchecked") public static void main(String[] args) throws CloneNotSupportedException { Ranks<Integer> matrix = new Matrix<>(new Integer[][]{{1,2,3},{3,4,5}}); System.out.println(matrix); System.out.println(Arrays.toString(matrix.shape())); System.out.println(matrix.rowNum()); System.out.println(matrix.colNum()); System.out.println(matrix.size()); System.out.println(matrix.get(4)); System.out.println(matrix.get(new int[]{0,0})); System.out.println(matrix.rowVector(1)); System.out.println(matrix.colVector(2)); Ranks<Integer> matrix1 = new Matrix<>(new Integer[][]{{7,8,9},{10,20,30}}); Addr<Integer> addr = new Addr<Integer>() { @Override public Integer zero() { return 0; } @Override public Integer add(Integer lhs, Integer rhs) { return lhs + rhs; } @Override public Integer sub(Integer lhs, Integer rhs) { return lhs - rhs; } @Override public Integer mul(Integer k, Integer hs) { return k * hs; } @Override public Integer square(Integer hs) { return hs * hs; } @Override public double prescription(Integer hs) { return Math.sqrt(hs); } @Override public double div(Integer hs, double nhs) { return hs / nhs; } @Override public double acos(double hs) { return Math.acos(hs); } @Override public double toDegree(double hs) { return Math.toDegrees(hs); } }; Ranks<Integer> copy = ((Matrix)matrix).clone(); Ranks<Integer> copy1 = ((Matrix)matrix).clone(); System.out.println(matrix.add(matrix1,addr)); System.out.println(matrix1.sub(copy,addr)); System.out.println(copy.div(2,addr)); System.out.println(copy.pos(addr,1)); System.out.println(copy.neg(addr,-1)); System.out.println(copy1.mul(3,addr)); System.out.println(Matrix.zero(2,3,addr)); } }
运行结果
Matrix{values=[[1, 2, 3],[3, 4, 5]]}
[2, 3]
2
3
6
4
1
Vector{values=[3, 4, 5]}
Vector{values=[3, 5]}
Matrix{values=[[8, 10, 12],[13, 24, 35]]}
Matrix{values=[[6, 6, 6],[7, 16, 25]]}
Matrix{values=[[0.5, 1.0, 1.5],[1.5, 2.0, 2.5]]}
Matrix{values=[[1, 2, 3],[3, 4, 5]]}
Matrix{values=[[-1, -2, -3],[-3, -4, -5]]}
Matrix{values=[[3, 6, 9],[9, 12, 15]]}
Matrix{values=[[0, 0, 0],[0, 0, 0]]}
Python代码
from .Vector import Vector class Matrix: def __init__(self, list2d): self._values = [row[:] for row in list2d] @classmethod def zero(cls, r, c): # 返回一个r行c列的零矩阵 return cls([[0] * c for _ in range(r)]) def __add__(self, another): # 返回两个矩阵的加法结果 assert self.shape() == another.shape(), \ "加法错误，两个矩阵的形状必须相同" return Matrix([[a + b for a, b in zip(self.row_vector(i), another.row_vector(i))] for i in range(self.row_num())]) def __sub__(self, another): # 返回两个矩阵的减法结果 assert self.shape() == another.shape(), \ "减法错误，两个矩阵的形状必须相同" return Matrix([[a - b for a, b in zip(self.row_vector(i), another.row_vector(i))] for i in range(self.row_num())]) def __mul__(self, k): # 返回矩阵的数量乘结果: self * k return Matrix([[e * k for e in self.row_vector(i)] for i in range(self.row_num())]) def __rmul__(self, k): # 返回矩阵的数量乘结果: k * self return self * k def __truediv__(self, k): # 返回数量除法的结果矩阵: self / k return (1 / k) * self def __pos__(self): # 返回矩阵取正的结果 return 1 * self def __neg__(self): # 返回矩阵取负的结果 return -1 * self def row_vector(self, index): # 返回矩阵的第index个行向量 return Vector(self._values[index]) def col_vector(self, index): # 返回矩阵的第index个列向量 return Vector([row[index] for row in self._values]) def __getitem__(self, pos): # 返回矩阵pos位置的元素 r, c = pos return self._values[r][c] def size(self): # 返回矩阵中元素的个数 r, c = self.shape() return r * c def row_num(self): # 返回矩阵的行数 return self.shape()[0] __len__ = row_num def col_num(self): # 返回矩阵的列数 return self.shape()[1] def shape(self): # 返回矩阵的形状：(行数，列数) return len(self._values), len(self._values[0]) def __repr__(self): return "Matrix({})".format(self._values) __str__ = __repr__
from playLA.Matrix import Matrix if __name__ == "__main__": matrix = Matrix([[1, 2, 3], [3, 4, 5]]) print(matrix) print("matrix.shape = {}".format(matrix.shape())) print("matrix.size = {}".format(matrix.size())) print("len(matrix) = {}".format(len(matrix))) print("matrix[0][0] = {}".format(matrix[0, 0])) print(matrix.row_vector(1)) print(matrix.col_vector(0)) matrix2 = Matrix([[7, 8, 9], [10, 20, 30]]) print("add: {}".format(matrix + matrix2)) print("subtract: {}".format(matrix2 - matrix)) print("scalar-mul: {}".format(matrix * 3)) print("scalar-mul: {}".format(3 * matrix)) print("zero_2_3: {}".format(Matrix.zero(2, 3)))
运行结果
Matrix([[1, 2, 3], [3, 4, 5]])
matrix.shape = (2, 3)
matrix.size = 6
len(matrix) = 2
matrix[0][0] = 1
(3, 4, 5)
(1, 3)
add: Matrix([[8, 10, 12], [13, 24, 35]])
subtract: Matrix([[6, 6, 6], [7, 16, 25]])
scalar-mul: Matrix([[3, 6, 9], [9, 12, 15]])
scalar-mul: Matrix([[3, 6, 9], [9, 12, 15]])
zero_2_3: Matrix([[0, 0, 0], [0, 0, 0]])
矩阵和向量的乘法
矩阵可以表示一个系统
比如在经济系统中，对IT、电子、矿产、房产的投入分别为
、
、
、
它们满足这样一个方程组
这里面每一个方程的意思为，如果我要发展一个行业，首先需要基础投资多少钱，比如要发展IT业，首先要基础投资100万亿，但是电子业对IT业有影响，为了发展IT业，需要对电子业多投资20%，需要对矿产业多投资10%，需要对房产业多投资50%。其他的方程的意思类同。这里是一个举例，不代表实际意义。
上面的方程组可以转化成下面这种形式
类似于上面这种方程组，我们称为线性方程组，线性方程组在各个领域，有着重要的作用。而这样的系统，在线性代数中，称为线性系统。
在上面的线形方程组中，我们将每一个未知数的系数提取出来，形成一个矩阵，就有了如下的矩阵
，而右侧部分的数提取出来，就形成了一个列向量
我们将该线形方程组的所有的未知数也全部提取出来，形成一个列向量
那么这个线形方程组，就等价于下面这种形式，系数矩阵乘以未知数列向量等于常数列向量
上述这个式子可以表示为
，其中A为矩阵，
是一个向量，
是一个向量。
因为该式子也之前的线形方程组等价，所以，该等式的左边可以用这样一个向量来表示
我们应该注意，这是一个向量，而不是一个矩阵。则我们可以得到如下的形式
根据这个等式，我们可以看出，矩阵与向量相乘
则是矩阵中的每一个行向量与向量进行点乘得到一个新的向量。我们之前知道，一个向量与一个向量进行点乘，得到的是一个标量，那么多个向量与一个向量进行点乘，就得到多个标量，再将这多个标量就可以组成一个向量。则我们得到一个法则
矩阵与向量相乘，就是该矩阵的行向量与该向量的点乘，得到一个新的向量。
转换成具体的公式，则为
我们知道向量与向量相乘，就是各个维度的分量值相乘再相加，但前提条件是这两个向量的长度必须向等。换成矩阵来说，就是矩阵的列数必须与向量的长度相等。即矩阵A的列数必须和向量u的元素个数一致。但矩阵A的行数没有限制。
换一种形式即为
这里需要说明的是，这里是从矩阵的行视角来看待矩阵的乘法的，后续会从矩阵的列视角来看待这一问题。
从来看，实际上就是矩阵A将向量x转换成了向量b,可以把矩阵理解成向量的函数。函数f(x)=b就是将x转换成b，而这里是向量x转换成了向量b。
首先，我们来扩展Ranks接口，新增一个矩阵与向量的点乘
public interface Ranks<T extends Number> { /** * 矩阵的行列数 * @return */ int[] shape(); /** * 矩阵的行数 * @return */ int rowNum(); /** * 矩阵的列数 * @return */ int colNum(); /** * 矩阵的元素个数 * @return */ int size(); /** * 获取矩阵中的一个元素 * @param pos * @return */ T get(int pos); /** * 获取矩阵中的一个元素 * @param pos * @return */ T get(int[] pos); /** * 获取一个行向量 * @param index * @return */ Attribute<T> rowVector(int index); /** * 获取一个列向量 * @param index * @return */ Attribute<T> colVector(int index); /** * 矩阵相加 * @param another * @return */ Ranks<T> add(Ranks<T> another,Addr<T> addr); /** * 矩阵相减 * @param another * @param addr * @return */ Ranks<T> sub(Ranks<T> another,Addr<T> addr); /** * 矩阵乘以一个数 * @param k * @param addr * @return */ Ranks<T> mul(T k,Addr<T> addr); /** * 矩阵除以一个数 * @param k * @param addr * @return */ Ranks<Double> div(double k,Addr<T> addr); /** * 矩阵取正的结果 * @return */ Ranks<T> pos(Addr<T> addr,T one); /** * 矩阵取负的结果 * @return */ Ranks<T> neg(Addr<T> addr,T negOne); /** * 矩阵与向量的点乘 * @param another * @param addr * @return */ Attribute<T> dot(Attribute<T> another,Addr<T> addr); }
矩阵实现类
public class Matrix<T extends Number> implements Ranks<T>,Cloneable { private T[][] values; public Matrix(T[][] values) { this.values = values; } @Override public String toString() { StringBuilder builder = new StringBuilder(); for (int i = 0; i < values.length - 1; i++) { builder.append(Arrays.toString(values[i]) + ","); } builder.append(Arrays.toString(values[values.length - 1])); return "Matrix{" + "values=" + "[" + builder.toString() + "]" + '}'; } /** * 创建一个零矩阵 * @param row * @param col * @param addr * @param <T> * @return */ @SuppressWarnings("unchecked") public static <T extends Number> Ranks<T> zero(int row,int col,Addr<T> addr) { T[][] values = (T[][])new Number[row][col]; for (int i = 0; i < row; i++) { Arrays.fill(values[i],addr.zero()); } return new Matrix<>(values); } @Override public int[] shape() { int[] res = new int[2]; res[0] = values.length; res[1] = values[0].length; return res; } @Override public int rowNum() { return shape()[0]; } @Override public int colNum() { return shape()[1]; } @Override public int size() { int row = rowNum(); int col = colNum(); return row * col; } @Override public T get(int pos) { if (pos >= size() || pos < 0) { throw new IllegalArgumentException("位置超出索引范围"); } int row = pos / colNum(); int col = pos % colNum(); return values[row][col]; } @Override public T get(int[] pos) { if (pos.length != 2) { throw new IllegalArgumentException("位置长度超出范围"); } if (pos[0] >= rowNum() || pos[0] < 0 || pos[1] >= colNum() || pos[1] < 0) { throw new IllegalArgumentException("位置超出索引范围"); } return values[pos[0]][pos[1]]; } @Override public Attribute<T> rowVector(int index) { if (index >= rowNum() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } return new Vector<>(values[index]); } @Override @SuppressWarnings("unchecked") public Attribute<T> colVector(int index) { if (index >= colNum() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } T[] values = (T[])new Number[rowNum()]; for (int i = 0; i < rowNum(); i++) { values[i] = this.values[i][index]; } return new Vector<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> add(Ranks<T> another,Addr<T> addr) { if (!Arrays.equals(this.shape(),another.shape())) { throw new IllegalArgumentException("加法错误，两个矩阵的形状必须相同"); } T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.add(this.get(new int[]{i,j}),another.get(new int[]{i,j})); } } this.values = values; return this; } @Override @SuppressWarnings("unchecked") public Ranks<T> sub(Ranks<T> another, Addr<T> addr) { if (!Arrays.equals(this.shape(),another.shape())) { throw new IllegalArgumentException("减法错误，两个矩阵的形状必须相同"); } T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.sub(this.get(new int[]{i,j}),another.get(new int[]{i,j})); } } this.values = values; return this; } @Override @SuppressWarnings("unchecked") public Ranks<T> mul(T k, Addr<T> addr) { T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.mul(k,this.get(new int[]{i,j})); } } this.values = values; return this; } @Override public Ranks<Double> div(double k, Addr<T> addr) { Double[][] values = new Double[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.div(this.get(new int[]{i,j}),k); } } return new Matrix<>(values); } @Override public Ranks<T> pos(Addr<T> addr,T one) { return mul(one,addr); } @Override public Ranks<T> neg(Addr<T> addr,T negOne) { return mul(negOne,addr); } @Override @SuppressWarnings("unchecked") public Attribute<T> dot(Attribute<T> another, Addr<T> addr) { if (this.colNum() != another.len()) { throw new IllegalArgumentException("矩阵列数需要与向量的长度相等"); } T[] values = (T[])new Number[this.rowNum()]; for (int i = 0; i < this.rowNum(); i++) { values[i] = this.rowVector(i).dot(another,addr); } return new Vector<>(values); } @Override public Matrix clone() throws CloneNotSupportedException { return (Matrix)super.clone(); } @SuppressWarnings("unchecked") public static void main(String[] args) throws CloneNotSupportedException { Ranks<Integer> matrix = new Matrix<>(new Integer[][]{{1,2,3},{3,4,5}}); System.out.println(matrix); System.out.println(Arrays.toString(matrix.shape())); System.out.println(matrix.rowNum()); System.out.println(matrix.colNum()); System.out.println(matrix.size()); System.out.println(matrix.get(4)); System.out.println(matrix.get(new int[]{0,0})); System.out.println(matrix.rowVector(1)); System.out.println(matrix.colVector(2)); Ranks<Integer> matrix1 = new Matrix<>(new Integer[][]{{7,8,9},{10,20,30}}); Addr<Integer> addr = new Addr<Integer>() { @Override public Integer zero() { return 0; } @Override public Integer add(Integer lhs, Integer rhs) { return lhs + rhs; } @Override public Integer sub(Integer lhs, Integer rhs) { return lhs - rhs; } @Override public Integer mul(Integer k, Integer hs) { return k * hs; } @Override public Integer square(Integer hs) { return hs * hs; } @Override public double prescription(Integer hs) { return Math.sqrt(hs); } @Override public double div(Integer hs, double nhs) { return hs / nhs; } @Override public double acos(double hs) { return Math.acos(hs); } @Override public double toDegree(double hs) { return Math.toDegrees(hs); } }; Ranks<Integer> copy = ((Matrix)matrix).clone(); Attribute<Integer> vector = new Vector<>(new Integer[]{23,65,79}); //[(1,2,3),(3,4,5)]*(23,65,79) System.out.println(matrix.dot(vector,addr)); Ranks<Integer> copy1 = ((Matrix)matrix).clone(); System.out.println(matrix.add(matrix1,addr)); System.out.println(matrix1.sub(copy,addr)); System.out.println(copy.div(2,addr)); System.out.println(copy.pos(addr,1)); System.out.println(copy.neg(addr,-1)); System.out.println(copy1.mul(3,addr)); System.out.println(Matrix.zero(2,3,addr)); } }
运行结果
Matrix{values=[[1, 2, 3],[3, 4, 5]]}
[2, 3]
2
3
6
4
1
Vector{values=[3, 4, 5]}
Vector{values=[3, 5]}
Vector{values=[390, 724]}
Matrix{values=[[8, 10, 12],[13, 24, 35]]}
Matrix{values=[[6, 6, 6],[7, 16, 25]]}
Matrix{values=[[0.5, 1.0, 1.5],[1.5, 2.0, 2.5]]}
Matrix{values=[[1, 2, 3],[3, 4, 5]]}
Matrix{values=[[-1, -2, -3],[-3, -4, -5]]}
Matrix{values=[[3, 6, 9],[9, 12, 15]]}
Matrix{values=[[0, 0, 0],[0, 0, 0]]}
Python代码
from .Vector import Vector class Matrix: def __init__(self, list2d): self._values = [row[:] for row in list2d] @classmethod def zero(cls, r, c): # 返回一个r行c列的零矩阵 return cls([[0] * c for _ in range(r)]) def __add__(self, another): # 返回两个矩阵的加法结果 assert self.shape() == another.shape(), \ "加法错误，两个矩阵的形状必须相同" return Matrix([[a + b for a, b in zip(self.row_vector(i), another.row_vector(i))] for i in range(self.row_num())]) def __sub__(self, another): # 返回两个矩阵的减法结果 assert self.shape() == another.shape(), \ "减法错误，两个矩阵的形状必须相同" return Matrix([[a - b for a, b in zip(self.row_vector(i), another.row_vector(i))] for i in range(self.row_num())]) def dot(self, another): if isinstance(another, Vector): # 矩阵和向量的乘法 assert self.col_num() == len(another), \ "矩阵列数需要与向量的长度相等" return Vector([self.row_vector(i).dot(another) for i in range(self.row_num())]) def __mul__(self, k): # 返回矩阵的数量乘结果: self * k return Matrix([[e * k for e in self.row_vector(i)] for i in range(self.row_num())]) def __rmul__(self, k): # 返回矩阵的数量乘结果: k * self return self * k def __truediv__(self, k): # 返回数量除法的结果矩阵: self / k return (1 / k) * self def __pos__(self): # 返回矩阵取正的结果 return 1 * self def __neg__(self): # 返回矩阵取负的结果 return -1 * self def row_vector(self, index): # 返回矩阵的第index个行向量 return Vector(self._values[index]) def col_vector(self, index): # 返回矩阵的第index个列向量 return Vector([row[index] for row in self._values]) def __getitem__(self, pos): # 返回矩阵pos位置的元素 r, c = pos return self._values[r][c] def size(self): # 返回矩阵中元素的个数 r, c = self.shape() return r * c def row_num(self): # 返回矩阵的行数 return self.shape()[0] __len__ = row_num def col_num(self): # 返回矩阵的列数 return self.shape()[1] def shape(self): # 返回矩阵的形状：(行数，列数) return len(self._values), len(self._values[0]) def __repr__(self): return "Matrix({})".format(self._values) __str__ = __repr__
from playLA.Matrix import Matrix from playLA.Vector import Vector if __name__ == "__main__": matrix = Matrix([[1, 2, 3], [3, 4, 5]]) print(matrix) print("matrix.shape = {}".format(matrix.shape())) print("matrix.size = {}".format(matrix.size())) print("len(matrix) = {}".format(len(matrix))) print("matrix[0][0] = {}".format(matrix[0, 0])) print(matrix.row_vector(1)) print(matrix.col_vector(0)) matrix2 = Matrix([[7, 8, 9], [10, 20, 30]]) print("add: {}".format(matrix + matrix2)) print("subtract: {}".format(matrix2 - matrix)) print("scalar-mul: {}".format(matrix * 3)) print("scalar-mul: {}".format(3 * matrix)) print("zero_2_3: {}".format(Matrix.zero(2, 3))) vec = Vector([23, 65, 79]) print("matrix.dot(vec) = {}".format(matrix.dot(vec)))
运行结果
Matrix([[1, 2, 3], [3, 4, 5]])
matrix.shape = (2, 3)
matrix.size = 6
len(matrix) = 2
matrix[0][0] = 1
(3, 4, 5)
(1, 3)
add: Matrix([[8, 10, 12], [13, 24, 35]])
subtract: Matrix([[6, 6, 6], [7, 16, 25]])
scalar-mul: Matrix([[3, 6, 9], [9, 12, 15]])
scalar-mul: Matrix([[3, 6, 9], [9, 12, 15]])
zero_2_3: Matrix([[0, 0, 0], [0, 0, 0]])
matrix.dot(vec) = (390, 724)
矩阵和矩阵的乘法
我们之前在说矩阵的数量乘法的时候，说到
则它乘以2之后为
但如果让每个点的横坐标扩大1.5倍，纵坐标扩大2倍
对于P中的每一个点其实是一个向量，只不过我们以列向量的形式来表现它
，其中的x可能为0、4、5；y可能为0、0、3
如此我们需要将其转换为
，由于矩阵是向量的一个函数，相当于我们要找到一个矩阵T满足
首先要满足矩阵与向量相乘，则T的列数一定与向量的长度相同，则T的列数一定是2。由于结果向量只有两行，所以T的行数也只有2.
由此我们可以得到下面这个式子，来满足矩阵乘以一个向量
因为结果为，由此我们可以很容易得出
最终我们可以得出
因为在矩阵P中有3个点，我们将每一个点都转成一个列向量，从而形成一个新的矩阵
这样一来，我们就可以让矩阵T与矩阵P相乘
来批量的把P中的点坐标变成另外一组点坐标
我们可以把P中的每一列都当成一个单独的向量，则对应的结果为
对应二维平面坐标中就是这个样子
由这个式子
我们可以看出矩阵和矩阵相乘，满足下面这个样子
矩阵和矩阵的乘法其实就是矩阵和向量乘法的批处理，由于矩阵和向量的乘法需要满足矩阵的列数必须和向量的长度相等，放在这里就是矩阵A的列数必须和矩阵B的行数相等。我们可以用这种形式来表示
其中的
到
就是一个一个的列向量，结果也是由矩阵A与这一个一个列向量点乘后的新的列向量所形成的新的矩阵。
最终得到的公式就是
如果A是m * k的矩阵，B是k * n的矩阵，则结果矩阵为m * n的矩阵。
矩阵乘法不遵守交换律 
，因为它们很有可能根本不能相乘，即便可以相乘，结果也不一样。
扩展Ranks接口，新增一个矩阵乘以矩阵的方法
public interface Ranks<T extends Number> { /** * 矩阵的行列数 * @return */ int[] shape(); /** * 矩阵的行数 * @return */ int rowNum(); /** * 矩阵的列数 * @return */ int colNum(); /** * 矩阵的元素个数 * @return */ int size(); /** * 获取矩阵中的一个元素 * @param pos * @return */ T get(int pos); /** * 获取矩阵中的一个元素 * @param pos * @return */ T get(int[] pos); /** * 获取一个行向量 * @param index * @return */ Attribute<T> rowVector(int index); /** * 获取一个列向量 * @param index * @return */ Attribute<T> colVector(int index); /** * 矩阵相加 * @param another * @return */ Ranks<T> add(Ranks<T> another,Addr<T> addr); /** * 矩阵相减 * @param another * @param addr * @return */ Ranks<T> sub(Ranks<T> another,Addr<T> addr); /** * 矩阵乘以一个数 * @param k * @param addr * @return */ Ranks<T> mul(T k,Addr<T> addr); /** * 矩阵除以一个数 * @param k * @param addr * @return */ Ranks<Double> div(double k,Addr<T> addr); /** * 矩阵取正的结果 * @return */ Ranks<T> pos(Addr<T> addr,T one); /** * 矩阵取负的结果 * @return */ Ranks<T> neg(Addr<T> addr,T negOne); /** * 矩阵与向量的点乘 * @param another * @param addr * @return */ Attribute<T> dot(Attribute<T> another,Addr<T> addr); /** * 矩阵与矩阵的点乘 * @param another * @param addr * @return */ Ranks<T> dot(Ranks<T> another,Addr<T> addr); }
矩阵实现类
public class Matrix<T extends Number> implements Ranks<T>,Cloneable { private T[][] values; public Matrix(T[][] values) { this.values = values; } @Override public String toString() { StringBuilder builder = new StringBuilder(); for (int i = 0; i < values.length - 1; i++) { builder.append(Arrays.toString(values[i]) + ","); } builder.append(Arrays.toString(values[values.length - 1])); return "Matrix{" + "values=" + "[" + builder.toString() + "]" + '}'; } /** * 创建一个零矩阵 * @param row * @param col * @param addr * @param <T> * @return */ @SuppressWarnings("unchecked") public static <T extends Number> Ranks<T> zero(int row,int col,Addr<T> addr) { T[][] values = (T[][])new Number[row][col]; for (int i = 0; i < row; i++) { Arrays.fill(values[i],addr.zero()); } return new Matrix<>(values); } @Override public int[] shape() { int[] res = new int[2]; res[0] = values.length; res[1] = values[0].length; return res; } @Override public int rowNum() { return shape()[0]; } @Override public int colNum() { return shape()[1]; } @Override public int size() { int row = rowNum(); int col = colNum(); return row * col; } @Override public T get(int pos) { if (pos >= size() || pos < 0) { throw new IllegalArgumentException("位置超出索引范围"); } int row = pos / colNum(); int col = pos % colNum(); return values[row][col]; } @Override public T get(int[] pos) { if (pos.length != 2) { throw new IllegalArgumentException("位置长度超出范围"); } if (pos[0] >= rowNum() || pos[0] < 0 || pos[1] >= colNum() || pos[1] < 0) { throw new IllegalArgumentException("位置超出索引范围"); } return values[pos[0]][pos[1]]; } @Override public Attribute<T> rowVector(int index) { if (index >= rowNum() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } return new Vector<>(values[index]); } @Override @SuppressWarnings("unchecked") public Attribute<T> colVector(int index) { if (index >= colNum() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } T[] values = (T[])new Number[rowNum()]; for (int i = 0; i < rowNum(); i++) { values[i] = this.values[i][index]; } return new Vector<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> add(Ranks<T> another,Addr<T> addr) { if (!Arrays.equals(this.shape(),another.shape())) { throw new IllegalArgumentException("加法错误，两个矩阵的形状必须相同"); } T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.add(this.get(new int[]{i,j}),another.get(new int[]{i,j})); } } this.values = values; return this; } @Override @SuppressWarnings("unchecked") public Ranks<T> sub(Ranks<T> another, Addr<T> addr) { if (!Arrays.equals(this.shape(),another.shape())) { throw new IllegalArgumentException("减法错误，两个矩阵的形状必须相同"); } T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.sub(this.get(new int[]{i,j}),another.get(new int[]{i,j})); } } this.values = values; return this; } @Override @SuppressWarnings("unchecked") public Ranks<T> mul(T k, Addr<T> addr) { T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.mul(k,this.get(new int[]{i,j})); } } this.values = values; return this; } @Override public Ranks<Double> div(double k, Addr<T> addr) { Double[][] values = new Double[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.div(this.get(new int[]{i,j}),k); } } return new Matrix<>(values); } @Override public Ranks<T> pos(Addr<T> addr,T one) { return mul(one,addr); } @Override public Ranks<T> neg(Addr<T> addr,T negOne) { return mul(negOne,addr); } @Override @SuppressWarnings("unchecked") public Attribute<T> dot(Attribute<T> another, Addr<T> addr) { if (this.colNum() != another.len()) { throw new IllegalArgumentException("矩阵列数需要与向量的长度相等"); } T[] values = (T[])new Number[this.rowNum()]; for (int i = 0; i < this.rowNum(); i++) { values[i] = this.rowVector(i).dot(another,addr); } return new Vector<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> dot(Ranks<T> another, Addr<T> addr) { if (this.colNum() != another.rowNum()) { throw new IllegalArgumentException("当前矩阵列数必须等于点乘矩阵行数"); } T[][] values = (T[][]) new Number[this.rowNum()][another.colNum()]; for (int i = 0; i < this.rowNum(); i++) { for (int j = 0; j < another.colNum(); j++) { values[i][j] = this.rowVector(i).dot(another.colVector(j),addr); } } return new Matrix<>(values); } @Override public Matrix clone() throws CloneNotSupportedException { return (Matrix)super.clone(); } @SuppressWarnings("unchecked") public static void main(String[] args) throws CloneNotSupportedException { Ranks<Integer> matrix = new Matrix<>(new Integer[][]{{1,2,3},{3,4,5}}); System.out.println(matrix); System.out.println(Arrays.toString(matrix.shape())); System.out.println(matrix.rowNum()); System.out.println(matrix.colNum()); System.out.println(matrix.size()); System.out.println(matrix.get(4)); System.out.println(matrix.get(new int[]{0,0})); System.out.println(matrix.rowVector(1)); System.out.println(matrix.colVector(2)); Ranks<Integer> matrix1 = new Matrix<>(new Integer[][]{{7,8,9},{10,20,30}}); Addr<Integer> addr = new Addr<Integer>() { @Override public Integer zero() { return 0; } @Override public Integer add(Integer lhs, Integer rhs) { return lhs + rhs; } @Override public Integer sub(Integer lhs, Integer rhs) { return lhs - rhs; } @Override public Integer mul(Integer k, Integer hs) { return k * hs; } @Override public Integer square(Integer hs) { return hs * hs; } @Override public double prescription(Integer hs) { return Math.sqrt(hs); } @Override public double div(Integer hs, double nhs) { return hs / nhs; } @Override public double acos(double hs) { return Math.acos(hs); } @Override public double toDegree(double hs) { return Math.toDegrees(hs); } }; Ranks<Integer> copy = ((Matrix)matrix).clone(); Attribute<Integer> vector = new Vector<>(new Integer[]{23,65,79}); Ranks<Integer> matrix2 = new Matrix<>(new Integer[][]{{7,9},{11,33},{42,97}}); //[(1,2,3),(3,4,5)]*(23,65,79) System.out.println(matrix.dot(vector,addr)); // [(1,2,3),(3,4,5)]*[(7,9),(11,33),(42,97)] System.out.println(matrix.dot(matrix2,addr)); Ranks<Integer> copy1 = ((Matrix)matrix).clone(); System.out.println(matrix.add(matrix1,addr)); System.out.println(matrix1.sub(copy,addr)); System.out.println(copy.div(2,addr)); System.out.println(copy.pos(addr,1)); System.out.println(copy.neg(addr,-1)); System.out.println(copy1.mul(3,addr)); System.out.println(Matrix.zero(2,3,addr)); } }
运行结果
Matrix{values=[[1, 2, 3],[3, 4, 5]]}
[2, 3]
2
3
6
4
1
Vector{values=[3, 4, 5]}
Vector{values=[3, 5]}
Vector{values=[390, 724]}
Matrix{values=[[155, 366],[275, 644]]}
Matrix{values=[[8, 10, 12],[13, 24, 35]]}
Matrix{values=[[6, 6, 6],[7, 16, 25]]}
Matrix{values=[[0.5, 1.0, 1.5],[1.5, 2.0, 2.5]]}
Matrix{values=[[1, 2, 3],[3, 4, 5]]}
Matrix{values=[[-1, -2, -3],[-3, -4, -5]]}
Matrix{values=[[3, 6, 9],[9, 12, 15]]}
Matrix{values=[[0, 0, 0],[0, 0, 0]]}
Python代码
from .Vector import Vector class Matrix: def __init__(self, list2d): self._values = [row[:] for row in list2d] @classmethod def zero(cls, r, c): # 返回一个r行c列的零矩阵 return cls([[0] * c for _ in range(r)]) def __add__(self, another): # 返回两个矩阵的加法结果 assert self.shape() == another.shape(), \ "加法错误，两个矩阵的形状必须相同" return Matrix([[a + b for a, b in zip(self.row_vector(i), another.row_vector(i))] for i in range(self.row_num())]) def __sub__(self, another): # 返回两个矩阵的减法结果 assert self.shape() == another.shape(), \ "减法错误，两个矩阵的形状必须相同" return Matrix([[a - b for a, b in zip(self.row_vector(i), another.row_vector(i))] for i in range(self.row_num())]) def dot(self, another): if isinstance(another, Vector): # 矩阵和向量的乘法 assert self.col_num() == len(another), \ "矩阵列数需要与向量的长度相等" return Vector([self.row_vector(i).dot(another) for i in range(self.row_num())]) if isinstance(another, Matrix): # 矩阵和矩阵的乘法 assert self.col_num() == another.row_num(), \ "当前矩阵列数必须等于点乘矩阵行数" return Matrix([[self.row_vector(i).dot(another.col_vector(j)) for j in range(another.col_num())] for i in range(self.row_num())]) def __mul__(self, k): # 返回矩阵的数量乘结果: self * k return Matrix([[e * k for e in self.row_vector(i)] for i in range(self.row_num())]) def __rmul__(self, k): # 返回矩阵的数量乘结果: k * self return self * k def __truediv__(self, k): # 返回数量除法的结果矩阵: self / k return (1 / k) * self def __pos__(self): # 返回矩阵取正的结果 return 1 * self def __neg__(self): # 返回矩阵取负的结果 return -1 * self def row_vector(self, index): # 返回矩阵的第index个行向量 return Vector(self._values[index]) def col_vector(self, index): # 返回矩阵的第index个列向量 return Vector([row[index] for row in self._values]) def __getitem__(self, pos): # 返回矩阵pos位置的元素 r, c = pos return self._values[r][c] def size(self): # 返回矩阵中元素的个数 r, c = self.shape() return r * c def row_num(self): # 返回矩阵的行数 return self.shape()[0] __len__ = row_num def col_num(self): # 返回矩阵的列数 return self.shape()[1] def shape(self): # 返回矩阵的形状：(行数，列数) return len(self._values), len(self._values[0]) def __repr__(self): return "Matrix({})".format(self._values) __str__ = __repr__
from playLA.Matrix import Matrix from playLA.Vector import Vector if __name__ == "__main__": matrix = Matrix([[1, 2, 3], [3, 4, 5]]) print(matrix) print("matrix.shape = {}".format(matrix.shape())) print("matrix.size = {}".format(matrix.size())) print("len(matrix) = {}".format(len(matrix))) print("matrix[0][0] = {}".format(matrix[0, 0])) print(matrix.row_vector(1)) print(matrix.col_vector(0)) matrix2 = Matrix([[7, 8, 9], [10, 20, 30]]) print("add: {}".format(matrix + matrix2)) print("subtract: {}".format(matrix2 - matrix)) print("scalar-mul: {}".format(matrix * 3)) print("scalar-mul: {}".format(3 * matrix)) print("zero_2_3: {}".format(Matrix.zero(2, 3))) vec = Vector([23, 65, 79]) print("matrix.dot(vec) = {}".format(matrix.dot(vec))) matrix3 = Matrix([[7, 9], [11, 33], [42, 97]]) print("matrix.dot(matrix3) = {}".format(matrix.dot(matrix3)))
运行结果
Matrix([[1, 2, 3], [3, 4, 5]])
matrix.shape = (2, 3)
matrix.size = 6
len(matrix) = 2
matrix[0][0] = 1
(3, 4, 5)
(1, 3)
add: Matrix([[8, 10, 12], [13, 24, 35]])
subtract: Matrix([[6, 6, 6], [7, 16, 25]])
scalar-mul: Matrix([[3, 6, 9], [9, 12, 15]])
scalar-mul: Matrix([[3, 6, 9], [9, 12, 15]])
zero_2_3: Matrix([[0, 0, 0], [0, 0, 0]])
matrix.dot(vec) = (390, 724)
matrix.dot(matrix3) = Matrix([[155, 366], [275, 644]])
矩阵乘法的性质
矩阵乘法不遵守交换律，但是矩阵乘法遵守下面的性质
乘法结合律 
乘法分配律 
因为矩阵乘法不遵守交换律，所以上面两条分配律是不等价的
对任意r*c的矩阵A,存在c*x的矩阵O,满足：
，这里需要注意的是，虽然
与
都是零矩阵，但它们的形状是不同的。
对任意r*c的矩阵A,存在x*r的矩阵O,满足:
矩阵的幂
矩阵A的k次方，也就是k个A相乘。由于矩阵相乘有条件限制，所以这里要做到k个A能够相乘，则A必须是一个方阵，只有方阵才可以进行矩阵的幂运算。
同样，由于矩阵的乘法不满足交换律，则
它最多可以扩展到
而中间的
不能合并成
矩阵的转置
之前我们说把变成这个样子的情况，称为矩阵的转置。因为在很多情况下，我们拿到的矩阵很有可能不是列矩阵，而是第一种情况的行矩阵，而为了方便我们做矩阵乘法，我们需要将第一种情况的行矩阵给变成列矩阵。
而我们对进行转置，可以写成
矩阵的转置：行变成列，列变成行。即为
，则
对于向量来说，有行向量
和列向量
，一般我们说向量如果不特殊声明的话，都是指一个列向量。由于横版印刷的原因，使用符号：
。
现在我们给Ranks接口增加一个矩阵转置的方法
public interface Ranks<T extends Number> { /** * 矩阵的行列数 * @return */ int[] shape(); /** * 矩阵的行数 * @return */ int rowNum(); /** * 矩阵的列数 * @return */ int colNum(); /** * 矩阵的元素个数 * @return */ int size(); /** * 获取矩阵中的一个元素 * @param pos * @return */ T get(int pos); /** * 获取矩阵中的一个元素 * @param pos * @return */ T get(int[] pos); /** * 获取一个行向量 * @param index * @return */ Attribute<T> rowVector(int index); /** * 获取一个列向量 * @param index * @return */ Attribute<T> colVector(int index); /** * 矩阵相加 * @param another * @return */ Ranks<T> add(Ranks<T> another,Addr<T> addr); /** * 矩阵相减 * @param another * @param addr * @return */ Ranks<T> sub(Ranks<T> another,Addr<T> addr); /** * 矩阵乘以一个数 * @param k * @param addr * @return */ Ranks<T> mul(T k,Addr<T> addr); /** * 矩阵除以一个数 * @param k * @param addr * @return */ Ranks<Double> div(double k,Addr<T> addr); /** * 矩阵取正的结果 * @return */ Ranks<T> pos(Addr<T> addr,T one); /** * 矩阵取负的结果 * @return */ Ranks<T> neg(Addr<T> addr,T negOne); /** * 矩阵与向量的点乘 * @param another * @param addr * @return */ Attribute<T> dot(Attribute<T> another,Addr<T> addr); /** * 矩阵与矩阵的点乘 * @param another * @param addr * @return */ Ranks<T> dot(Ranks<T> another,Addr<T> addr); /** * 转置矩阵 * @return */ Ranks<T> trans(); }
矩阵实现类
public class Matrix<T extends Number> implements Ranks<T>,Cloneable { private T[][] values; public Matrix(T[][] values) { this.values = values; } @Override public String toString() { StringBuilder builder = new StringBuilder(); for (int i = 0; i < values.length - 1; i++) { builder.append(Arrays.toString(values[i]) + ","); } builder.append(Arrays.toString(values[values.length - 1])); return "Matrix{" + "values=" + "[" + builder.toString() + "]" + '}'; } /** * 创建一个零矩阵 * @param row * @param col * @param addr * @param <T> * @return */ @SuppressWarnings("unchecked") public static <T extends Number> Ranks<T> zero(int row,int col,Addr<T> addr) { T[][] values = (T[][])new Number[row][col]; for (int i = 0; i < row; i++) { Arrays.fill(values[i],addr.zero()); } return new Matrix<>(values); } @Override public int[] shape() { int[] res = new int[2]; res[0] = values.length; res[1] = values[0].length; return res; } @Override public int rowNum() { return shape()[0]; } @Override public int colNum() { return shape()[1]; } @Override public int size() { int row = rowNum(); int col = colNum(); return row * col; } @Override public T get(int pos) { if (pos >= size() || pos < 0) { throw new IllegalArgumentException("位置超出索引范围"); } int row = pos / colNum(); int col = pos % colNum(); return values[row][col]; } @Override public T get(int[] pos) { if (pos.length != 2) { throw new IllegalArgumentException("位置长度超出范围"); } if (pos[0] >= rowNum() || pos[0] < 0 || pos[1] >= colNum() || pos[1] < 0) { throw new IllegalArgumentException("位置超出索引范围"); } return values[pos[0]][pos[1]]; } @Override public Attribute<T> rowVector(int index) { if (index >= rowNum() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } return new Vector<>(values[index]); } @Override @SuppressWarnings("unchecked") public Attribute<T> colVector(int index) { if (index >= colNum() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } T[] values = (T[])new Number[rowNum()]; for (int i = 0; i < rowNum(); i++) { values[i] = this.values[i][index]; } return new Vector<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> add(Ranks<T> another,Addr<T> addr) { if (!Arrays.equals(this.shape(),another.shape())) { throw new IllegalArgumentException("加法错误，两个矩阵的形状必须相同"); } T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.add(this.get(new int[]{i,j}),another.get(new int[]{i,j})); } } this.values = values; return this; } @Override @SuppressWarnings("unchecked") public Ranks<T> sub(Ranks<T> another, Addr<T> addr) { if (!Arrays.equals(this.shape(),another.shape())) { throw new IllegalArgumentException("减法错误，两个矩阵的形状必须相同"); } T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.sub(this.get(new int[]{i,j}),another.get(new int[]{i,j})); } } this.values = values; return this; } @Override @SuppressWarnings("unchecked") public Ranks<T> mul(T k, Addr<T> addr) { T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.mul(k,this.get(new int[]{i,j})); } } this.values = values; return this; } @Override public Ranks<Double> div(double k, Addr<T> addr) { Double[][] values = new Double[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.div(this.get(new int[]{i,j}),k); } } return new Matrix<>(values); } @Override public Ranks<T> pos(Addr<T> addr,T one) { return mul(one,addr); } @Override public Ranks<T> neg(Addr<T> addr,T negOne) { return mul(negOne,addr); } @Override @SuppressWarnings("unchecked") public Attribute<T> dot(Attribute<T> another, Addr<T> addr) { if (this.colNum() != another.len()) { throw new IllegalArgumentException("矩阵列数需要与向量的长度相等"); } T[] values = (T[])new Number[this.rowNum()]; for (int i = 0; i < this.rowNum(); i++) { values[i] = this.rowVector(i).dot(another,addr); } return new Vector<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> dot(Ranks<T> another, Addr<T> addr) { if (this.colNum() != another.rowNum()) { throw new IllegalArgumentException("当前矩阵列数必须等于点乘矩阵行数"); } T[][] values = (T[][]) new Number[this.rowNum()][another.colNum()]; for (int i = 0; i < this.rowNum(); i++) { for (int j = 0; j < another.colNum(); j++) { values[i][j] = this.rowVector(i).dot(another.colVector(j),addr); } } return new Matrix<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> trans() { T[][] values = (T[][])new Number[this.colNum()][this.rowNum()]; for (int i = 0; i < this.colNum(); i++) { for (int j = 0; j < this.rowNum(); j++) { values[i][j] = this.values[j][i]; } } return new Matrix<>(values); } @Override public Matrix clone() throws CloneNotSupportedException { return (Matrix)super.clone(); } @SuppressWarnings("unchecked") public static void main(String[] args) throws CloneNotSupportedException { Ranks<Integer> matrix = new Matrix<>(new Integer[][]{{1,2,3},{3,4,5}}); System.out.println(matrix); System.out.println(Arrays.toString(matrix.shape())); System.out.println(matrix.rowNum()); System.out.println(matrix.colNum()); System.out.println(matrix.size()); System.out.println(matrix.get(4)); System.out.println(matrix.get(new int[]{0,0})); System.out.println(matrix.rowVector(1)); System.out.println(matrix.colVector(2)); Ranks<Integer> matrix1 = new Matrix<>(new Integer[][]{{7,8,9},{10,20,30}}); Addr<Integer> addr = new Addr<Integer>() { @Override public Integer zero() { return 0; } @Override public Integer add(Integer lhs, Integer rhs) { return lhs + rhs; } @Override public Integer sub(Integer lhs, Integer rhs) { return lhs - rhs; } @Override public Integer mul(Integer k, Integer hs) { return k * hs; } @Override public Integer square(Integer hs) { return hs * hs; } @Override public double prescription(Integer hs) { return Math.sqrt(hs); } @Override public double div(Integer hs, double nhs) { return hs / nhs; } @Override public double acos(double hs) { return Math.acos(hs); } @Override public double toDegree(double hs) { return Math.toDegrees(hs); } }; System.out.println(matrix.trans()); Ranks<Integer> copy = ((Matrix)matrix).clone(); Attribute<Integer> vector = new Vector<>(new Integer[]{9,33,97}); Ranks<Integer> matrix2 = new Matrix<>(new Integer[][]{{7,9},{11,33},{42,97}}); //[(1,2,3),(3,4,5)]*(23,65,79) System.out.println(matrix.dot(vector,addr)); // [(1,2,3),(3,4,5)]*[(7,9),(11,33),(42,97)] System.out.println(matrix.dot(matrix2,addr)); Ranks<Integer> copy1 = ((Matrix)matrix).clone(); System.out.println(matrix.add(matrix1,addr)); System.out.println(matrix1.sub(copy,addr)); System.out.println(copy.div(2,addr)); System.out.println(copy.pos(addr,1)); System.out.println(copy.neg(addr,-1)); System.out.println(copy1.mul(3,addr)); System.out.println(Matrix.zero(2,3,addr)); } }
运行结果
Matrix{values=[[1, 2, 3],[3, 4, 5]]}
[2, 3]
2
3
6
4
1
Vector{values=[3, 4, 5]}
Vector{values=[3, 5]}
Matrix{values=[[1, 3],[2, 4],[3, 5]]}
Vector{values=[366, 644]}
Matrix{values=[[155, 366],[275, 644]]}
Matrix{values=[[8, 10, 12],[13, 24, 35]]}
Matrix{values=[[6, 6, 6],[7, 16, 25]]}
Matrix{values=[[0.5, 1.0, 1.5],[1.5, 2.0, 2.5]]}
Matrix{values=[[1, 2, 3],[3, 4, 5]]}
Matrix{values=[[-1, -2, -3],[-3, -4, -5]]}
Matrix{values=[[3, 6, 9],[9, 12, 15]]}
Matrix{values=[[0, 0, 0],[0, 0, 0]]}
Python代码
from .Vector import Vector class Matrix: def __init__(self, list2d): self._values = [row[:] for row in list2d] @classmethod def zero(cls, r, c): # 返回一个r行c列的零矩阵 return cls([[0] * c for _ in range(r)]) def trans(self): # 返回矩阵的转置矩阵 return Matrix([[e for e in self.col_vector(i)] for i in range(self.col_num())]) def __add__(self, another): # 返回两个矩阵的加法结果 assert self.shape() == another.shape(), \ "加法错误，两个矩阵的形状必须相同" return Matrix([[a + b for a, b in zip(self.row_vector(i), another.row_vector(i))] for i in range(self.row_num())]) def __sub__(self, another): # 返回两个矩阵的减法结果 assert self.shape() == another.shape(), \ "减法错误，两个矩阵的形状必须相同" return Matrix([[a - b for a, b in zip(self.row_vector(i), another.row_vector(i))] for i in range(self.row_num())]) def dot(self, another): if isinstance(another, Vector): # 矩阵和向量的乘法 assert self.col_num() == len(another), \ "矩阵列数需要与向量的长度相等" return Vector([self.row_vector(i).dot(another) for i in range(self.row_num())]) if isinstance(another, Matrix): # 矩阵和矩阵的乘法 assert self.col_num() == another.row_num(), \ "当前矩阵列数必须等于点乘矩阵行数" return Matrix([[self.row_vector(i).dot(another.col_vector(j)) for j in range(another.col_num())] for i in range(self.row_num())]) def __mul__(self, k): # 返回矩阵的数量乘结果: self * k return Matrix([[e * k for e in self.row_vector(i)] for i in range(self.row_num())]) def __rmul__(self, k): # 返回矩阵的数量乘结果: k * self return self * k def __truediv__(self, k): # 返回数量除法的结果矩阵: self / k return (1 / k) * self def __pos__(self): # 返回矩阵取正的结果 return 1 * self def __neg__(self): # 返回矩阵取负的结果 return -1 * self def row_vector(self, index): # 返回矩阵的第index个行向量 return Vector(self._values[index]) def col_vector(self, index): # 返回矩阵的第index个列向量 return Vector([row[index] for row in self._values]) def __getitem__(self, pos): # 返回矩阵pos位置的元素 r, c = pos return self._values[r][c] def size(self): # 返回矩阵中元素的个数 r, c = self.shape() return r * c def row_num(self): # 返回矩阵的行数 return self.shape()[0] __len__ = row_num def col_num(self): # 返回矩阵的列数 return self.shape()[1] def shape(self): # 返回矩阵的形状：(行数，列数) return len(self._values), len(self._values[0]) def __repr__(self): return "Matrix({})".format(self._values) __str__ = __repr__
from playLA.Matrix import Matrix from playLA.Vector import Vector if __name__ == "__main__": matrix = Matrix([[1, 2, 3], [3, 4, 5]]) print(matrix) print("matrix.shape = {}".format(matrix.shape())) print("matrix.size = {}".format(matrix.size())) print("len(matrix) = {}".format(len(matrix))) print("matrix[0][0] = {}".format(matrix[0, 0])) print(matrix.row_vector(1)) print(matrix.col_vector(0)) matrix2 = Matrix([[7, 8, 9], [10, 20, 30]]) print("add: {}".format(matrix + matrix2)) print("subtract: {}".format(matrix2 - matrix)) print("scalar-mul: {}".format(matrix * 3)) print("scalar-mul: {}".format(3 * matrix)) print("zero_2_3: {}".format(Matrix.zero(2, 3))) vec = Vector([23, 65, 79]) print("matrix.dot(vec) = {}".format(matrix.dot(vec))) matrix3 = Matrix([[7, 9], [11, 33], [42, 97]]) print("matrix.dot(matrix3) = {}".format(matrix.dot(matrix3))) print("matrix.trans = {}".format(matrix.trans()))
运行结果
Matrix([[1, 2, 3], [3, 4, 5]])
matrix.shape = (2, 3)
matrix.size = 6
len(matrix) = 2
matrix[0][0] = 1
(3, 4, 5)
(1, 3)
add: Matrix([[8, 10, 12], [13, 24, 35]])
subtract: Matrix([[6, 6, 6], [7, 16, 25]])
scalar-mul: Matrix([[3, 6, 9], [9, 12, 15]])
scalar-mul: Matrix([[3, 6, 9], [9, 12, 15]])
zero_2_3: Matrix([[0, 0, 0], [0, 0, 0]])
matrix.dot(vec) = (390, 724)
matrix.dot(matrix3) = Matrix([[155, 366], [275, 644]])
matrix.trans = Matrix([[1, 3], [2, 4], [3, 5]])
矩阵转置的性质
矩阵A的转置的转置就等于A本身 
两个矩阵相加后的转置等于分别对两个矩阵做转置后再相加 
矩阵数量乘后的转置等于数量乘以转置后的矩阵 
两个矩阵相乘后的转置等于第二个矩阵的转置乘以第一个矩阵的转置 
A是m*k的矩阵，B是k*n的矩阵，则A的转置就是k*m的矩阵，B的转置就是n*k的矩阵，则AB是m*n的矩阵，AB的转置就是n*m的矩阵。则(B的转置)(A的转置)是n*m的矩阵，满足上面的式子。
Numpy中矩阵的基本操作
import numpy as np if __name__ == "__main__": A = np.array([[1, 2], [3, 4]]) print(A) # 矩阵的形状 print(A.shape) # 矩阵的转置 print(A.T) # 获取矩阵的元素 print(A[1, 1]) # 矩阵的行向量 print(A[0]) print(A[1, :]) # 矩阵的列向量 print(A[:, 0]) B = np.array([[5, 6], [7, 8]]) # 矩阵加减法 print(A + B) print(A - B) # 矩阵的数量乘 print(10 * A) print(A * 10) # 矩阵的*号乘法，没有数学意义 print(A * B) # 矩阵与矩阵的点乘 print(A.dot(B)) vec = np.array([10, 100]) # 矩阵与向量的加法，没有数学意义 print(A + vec) print(A + 1) # 矩阵与向量的点乘 print(A.dot(vec))
运行结果
[[1 2]
[3 4]]
(2, 2)
[[1 3]
[2 4]]
4
[1 2]
[3 4]
[1 3]
[[ 6 8]
[10 12]]
[[-4 -4]
[-4 -4]]
[[10 20]
[30 40]]
[[10 20]
[30 40]]
[[ 5 12]
[21 32]]
[[19 22]
[43 50]]
[[ 11 102]
[ 13 104]]
[[2 3]
[4 5]]
[210 430]
图形变换矩阵(向量的函数)
之前我们知道了在直角坐标系中，让横坐标扩大1.5倍，纵坐标扩大2倍，则该变换矩阵是
则可以得出，让每个点到横坐标扩大a倍，纵坐标扩大b倍
因为矩阵与一个向量相乘，可得
现在我们来看一下更多的图形变换矩阵
- 让每个点关于x轴翻转
换句话说，即
，我们让
，代入进去，即
，这样我们就得到了
最终结果就是
- 让每个点关于y轴翻转
即
最终结果为
- 让每个点关于原点翻转(x轴，y轴均翻转)
即
最终结果为
我们也可以使用之前的结果来进行推导
沿x轴翻转的矩阵为
沿y轴翻转的矩阵为
则我们只需要把沿x轴翻转之后的点再沿y轴翻转就得到原点翻转的结果
- 沿x轴方向错切
即y轴不变，沿x轴进行变换，相当于
，最终得到
- 沿y方向的错切
即x轴不变，沿y轴进行变换，相当于
,最终得到
- 旋转
即坐标点沿着某个方向旋转了
角，这道题的关键在于需要知道(x,y)旋转之后的坐标
即我们需要知道
是多少？
由于一个向量旋转后，它的模是不变的(用d表示模)，则我们有
最终可以得到
根据两个角的差的余弦公式，可得
根据两个角的差的正弦公式，可得
这样，我们就知道了旋转后的坐标
最终可得
Python代码实现
我们先画一个图形
import matplotlib.pyplot as plt if __name__ == "__main__": points = [[0, 0], [0, 5], [3, 5], [3, 4], [1, 4], [1, 3], [2, 3], [2, 2], [1, 2], [1, 0]] x = [point[0] for point in points] y = [point[1] for point in points] plt.figure(figsize=(5, 5)) plt.xlim(-10, 10) plt.ylim(-10, 10) plt.plot(x, y) plt.show()
它就是这个样子的
现在我们要将该图形沿着x轴扩大2倍，沿着y轴扩大1.5倍
import matplotlib.pyplot as plt from playLA.Matrix import Matrix if __name__ == "__main__": points = [[0, 0], [0, 5], [3, 5], [3, 4], [1, 4], [1, 3], [2, 3], [2, 2], [1, 2], [1, 0]] x = [point[0] for point in points] y = [point[1] for point in points] plt.figure(figsize=(5, 5)) plt.xlim(-10,10) plt.ylim(-10,10) plt.plot(x,y) # plt.show() P = Matrix(points) # 缩放矩阵 T = Matrix([[2, 0], [0, 1.5]]) P2 = T.dot(P.trans()) plt.plot([P2.col_vector(i)[0] for i in range(P2.col_num())], [P2.col_vector(i)[1] for i in range(P2.col_num())]) plt.show()
得到的图形如下
现在我们让该图沿着x轴翻转
import matplotlib.pyplot as plt from playLA.Matrix import Matrix if __name__ == "__main__": points = [[0, 0], [0, 5], [3, 5], [3, 4], [1, 4], [1, 3], [2, 3], [2, 2], [1, 2], [1, 0]] x = [point[0] for point in points] y = [point[1] for point in points] plt.figure(figsize=(5, 5)) plt.xlim(-10,10) plt.ylim(-10,10) plt.plot(x,y) # plt.show() P = Matrix(points) # 缩放矩阵 # T = Matrix([[2, 0], [0, 1.5]]) # x轴翻转矩阵 T = Matrix([[1, 0], [0, -1]]) P2 = T.dot(P.trans()) plt.plot([P2.col_vector(i)[0] for i in range(P2.col_num())], [P2.col_vector(i)[1] for i in range(P2.col_num())]) plt.show()
得到的图形如下
现在我们让该图沿着y轴翻转
import matplotlib.pyplot as plt from playLA.Matrix import Matrix if __name__ == "__main__": points = [[0, 0], [0, 5], [3, 5], [3, 4], [1, 4], [1, 3], [2, 3], [2, 2], [1, 2], [1, 0]] x = [point[0] for point in points] y = [point[1] for point in points] plt.figure(figsize=(5, 5)) plt.xlim(-10,10) plt.ylim(-10,10) plt.plot(x,y) # plt.show() P = Matrix(points) # 缩放矩阵 # T = Matrix([[2, 0], [0, 1.5]]) # x轴翻转矩阵 # T = Matrix([[1, 0], [0, -1]]) # y轴翻转矩阵 T = Matrix([[-1, 0], [0, 1]]) P2 = T.dot(P.trans()) plt.plot([P2.col_vector(i)[0] for i in range(P2.col_num())], [P2.col_vector(i)[1] for i in range(P2.col_num())]) plt.show()
得到的图形如下
现在我们让该图沿着原点翻转
import matplotlib.pyplot as plt from playLA.Matrix import Matrix if __name__ == "__main__": points = [[0, 0], [0, 5], [3, 5], [3, 4], [1, 4], [1, 3], [2, 3], [2, 2], [1, 2], [1, 0]] x = [point[0] for point in points] y = [point[1] for point in points] plt.figure(figsize=(5, 5)) plt.xlim(-10,10) plt.ylim(-10,10) plt.plot(x,y) # plt.show() P = Matrix(points) # 缩放矩阵 # T = Matrix([[2, 0], [0, 1.5]]) # x轴翻转矩阵 # T = Matrix([[1, 0], [0, -1]]) # y轴翻转矩阵 # T = Matrix([[-1, 0], [0, 1]]) # 完全翻转矩阵 T = Matrix([[-1, 0], [0, -1]]) P2 = T.dot(P.trans()) plt.plot([P2.col_vector(i)[0] for i in range(P2.col_num())], [P2.col_vector(i)[1] for i in range(P2.col_num())]) plt.show()
得到的图形如下
现在我们让该图沿着x轴错切
import matplotlib.pyplot as plt from playLA.Matrix import Matrix if __name__ == "__main__": points = [[0, 0], [0, 5], [3, 5], [3, 4], [1, 4], [1, 3], [2, 3], [2, 2], [1, 2], [1, 0]] x = [point[0] for point in points] y = [point[1] for point in points] plt.figure(figsize=(5, 5)) plt.xlim(-10,10) plt.ylim(-10,10) plt.plot(x,y) # plt.show() P = Matrix(points) # 缩放矩阵 # T = Matrix([[2, 0], [0, 1.5]]) # x轴翻转矩阵 # T = Matrix([[1, 0], [0, -1]]) # y轴翻转矩阵 # T = Matrix([[-1, 0], [0, 1]]) # 完全翻转矩阵 # T = Matrix([[-1, 0], [0, -1]]) # x轴错切矩阵 T = Matrix([[1, 0.5], [0, 1]]) P2 = T.dot(P.trans()) plt.plot([P2.col_vector(i)[0] for i in range(P2.col_num())], [P2.col_vector(i)[1] for i in range(P2.col_num())]) plt.show()
得到的图形如下
现在我们让该图沿着y轴错切
import matplotlib.pyplot as plt from playLA.Matrix import Matrix if __name__ == "__main__": points = [[0, 0], [0, 5], [3, 5], [3, 4], [1, 4], [1, 3], [2, 3], [2, 2], [1, 2], [1, 0]] x = [point[0] for point in points] y = [point[1] for point in points] plt.figure(figsize=(5, 5)) plt.xlim(-10,10) plt.ylim(-10,10) plt.plot(x,y) # plt.show() P = Matrix(points) # 缩放矩阵 # T = Matrix([[2, 0], [0, 1.5]]) # x轴翻转矩阵 # T = Matrix([[1, 0], [0, -1]]) # y轴翻转矩阵 # T = Matrix([[-1, 0], [0, 1]]) # 完全翻转矩阵 # T = Matrix([[-1, 0], [0, -1]]) # x轴错切矩阵 # T = Matrix([[1, 0.5], [0, 1]]) # y轴错切矩阵 T = Matrix([[1, 0], [0.5, 1]]) P2 = T.dot(P.trans()) plt.plot([P2.col_vector(i)[0] for i in range(P2.col_num())], [P2.col_vector(i)[1] for i in range(P2.col_num())]) plt.show()
得到的图形如下
现在我们让该图旋转
import matplotlib.pyplot as plt from playLA.Matrix import Matrix import math if __name__ == "__main__": points = [[0, 0], [0, 5], [3, 5], [3, 4], [1, 4], [1, 3], [2, 3], [2, 2], [1, 2], [1, 0]] x = [point[0] for point in points] y = [point[1] for point in points] plt.figure(figsize=(5, 5)) plt.xlim(-10,10) plt.ylim(-10,10) plt.plot(x,y) # plt.show() P = Matrix(points) # 缩放矩阵 # T = Matrix([[2, 0], [0, 1.5]]) # x轴翻转矩阵 # T = Matrix([[1, 0], [0, -1]]) # y轴翻转矩阵 # T = Matrix([[-1, 0], [0, 1]]) # 完全翻转矩阵 # T = Matrix([[-1, 0], [0, -1]]) # x轴错切矩阵 # T = Matrix([[1, 0.5], [0, 1]]) # y轴错切矩阵 # T = Matrix([[1, 0], [0.5, 1]]) # 旋转操作 # 60度 theta = math.pi / 3 T = Matrix([[math.cos(theta), math.sin(theta)], [-math.sin(theta), math.cos(theta)]]) P2 = T.dot(P.trans()) plt.plot([P2.col_vector(i)[0] for i in range(P2.col_num())], [P2.col_vector(i)[1] for i in range(P2.col_num())]) plt.show()
得到的图形如下
单位矩阵
在我们的图形变换矩阵中，有一种特殊的变换矩阵，它不会改变其他图形的任何形状，即让每个点的横坐标保持1倍，纵坐标也保持1倍。
则该矩阵类似于这个样子
，当然这是在二维的平面直角坐标中
通常我们称这样的矩阵为单位矩阵
其中的2表示是二维的，而在三维空间中，单位矩阵则为
，同理在N维空间中单位矩阵就是
，又或者它可以记为
由以上可知，单位矩阵一定是方阵
单位矩阵的性质
则可得出
但是这里需要注意的是这个A不一定是一个方阵，如果A不为方阵，则上面的单位矩阵I的维度则不同，所以上面的等式并不是一个简单的乘法交换律，比如
如果要将式子左边的两个矩阵交换位置，则要满足矩阵乘法的条件(第一个矩阵的列数必须等于第二个矩阵的行数)，需要变成
所以这个单位矩阵需要变成一个三维的单位矩阵，但结果依然不变
我们需要在Addr接口中增加一个给出1的方法
public interface Addr<T extends Number> { T zero(); T one(); T add(T lhs,T rhs); T sub(T lhs,T rhs); T mul(T k,T hs); T square(T hs); double prescription(T hs); double div(T hs,double nhs); double acos(double hs); double toDegree(double hs); }
矩阵实现类，增加一个单位矩阵的静态方法
public class Matrix<T extends Number> implements Ranks<T>,Cloneable { private T[][] values; public Matrix(T[][] values) { this.values = values; } @Override public String toString() { StringBuilder builder = new StringBuilder(); for (int i = 0; i < values.length - 1; i++) { builder.append(Arrays.toString(values[i]) + ","); } builder.append(Arrays.toString(values[values.length - 1])); return "Matrix{" + "values=" + "[" + builder.toString() + "]" + '}'; } /** * 创建一个零矩阵 * @param row * @param col * @param addr * @param <T> * @return */ @SuppressWarnings("unchecked") public static <T extends Number> Ranks<T> zero(int row,int col,Addr<T> addr) { T[][] values = (T[][])new Number[row][col]; for (int i = 0; i < row; i++) { Arrays.fill(values[i],addr.zero()); } return new Matrix<>(values); } /** * 创建一个单位矩阵 * @param n * @param <T> * @return */ @SuppressWarnings("unchecked") public static <T extends Number> Ranks<T> identity(int n,Addr<T> addr) { T[][] values = (T[][])new Number[n][n]; for (int i = 0; i < n; i++) { Arrays.fill(values[i],addr.zero()); } for (int i = 0; i < n; i++) { values[i][i] = addr.one(); } return new Matrix<>(values); } @Override public int[] shape() { int[] res = new int[2]; res[0] = values.length; res[1] = values[0].length; return res; } @Override public int rowNum() { return shape()[0]; } @Override public int colNum() { return shape()[1]; } @Override public int size() { int row = rowNum(); int col = colNum(); return row * col; } @Override public T get(int pos) { if (pos >= size() || pos < 0) { throw new IllegalArgumentException("位置超出索引范围"); } int row = pos / colNum(); int col = pos % colNum(); return values[row][col]; } @Override public T get(int[] pos) { if (pos.length != 2) { throw new IllegalArgumentException("位置长度超出范围"); } if (pos[0] >= rowNum() || pos[0] < 0 || pos[1] >= colNum() || pos[1] < 0) { throw new IllegalArgumentException("位置超出索引范围"); } return values[pos[0]][pos[1]]; } @Override public Attribute<T> rowVector(int index) { if (index >= rowNum() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } return new Vector<>(values[index]); } @Override @SuppressWarnings("unchecked") public Attribute<T> colVector(int index) { if (index >= colNum() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } T[] values = (T[])new Number[rowNum()]; for (int i = 0; i < rowNum(); i++) { values[i] = this.values[i][index]; } return new Vector<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> add(Ranks<T> another,Addr<T> addr) { if (!Arrays.equals(this.shape(),another.shape())) { throw new IllegalArgumentException("加法错误，两个矩阵的形状必须相同"); } T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.add(this.get(new int[]{i,j}),another.get(new int[]{i,j})); } } this.values = values; return this; } @Override @SuppressWarnings("unchecked") public Ranks<T> sub(Ranks<T> another, Addr<T> addr) { if (!Arrays.equals(this.shape(),another.shape())) { throw new IllegalArgumentException("减法错误，两个矩阵的形状必须相同"); } T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.sub(this.get(new int[]{i,j}),another.get(new int[]{i,j})); } } this.values = values; return this; } @Override @SuppressWarnings("unchecked") public Ranks<T> mul(T k, Addr<T> addr) { T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.mul(k,this.get(new int[]{i,j})); } } this.values = values; return this; } @Override public Ranks<Double> div(double k, Addr<T> addr) { Double[][] values = new Double[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.div(this.get(new int[]{i,j}),k); } } return new Matrix<>(values); } @Override public Ranks<T> pos(Addr<T> addr,T one) { return mul(one,addr); } @Override public Ranks<T> neg(Addr<T> addr,T negOne) { return mul(negOne,addr); } @Override @SuppressWarnings("unchecked") public Attribute<T> dot(Attribute<T> another, Addr<T> addr) { if (this.colNum() != another.len()) { throw new IllegalArgumentException("矩阵列数需要与向量的长度相等"); } T[] values = (T[])new Number[this.rowNum()]; for (int i = 0; i < this.rowNum(); i++) { values[i] = this.rowVector(i).dot(another,addr); } return new Vector<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> dot(Ranks<T> another, Addr<T> addr) { if (this.colNum() != another.rowNum()) { throw new IllegalArgumentException("当前矩阵列数必须等于点乘矩阵行数"); } T[][] values = (T[][]) new Number[this.rowNum()][another.colNum()]; for (int i = 0; i < this.rowNum(); i++) { for (int j = 0; j < another.colNum(); j++) { values[i][j] = this.rowVector(i).dot(another.colVector(j),addr); } } return new Matrix<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> trans() { T[][] values = (T[][])new Number[this.colNum()][this.rowNum()]; for (int i = 0; i < this.colNum(); i++) { for (int j = 0; j < this.rowNum(); j++) { values[i][j] = this.values[j][i]; } } return new Matrix<>(values); } @Override public Matrix clone() throws CloneNotSupportedException { return (Matrix)super.clone(); } @SuppressWarnings("unchecked") public static void main(String[] args) throws CloneNotSupportedException { Ranks<Integer> matrix = new Matrix<>(new Integer[][]{{1,2,3},{3,4,5}}); System.out.println(matrix); System.out.println(Arrays.toString(matrix.shape())); System.out.println(matrix.rowNum()); System.out.println(matrix.colNum()); System.out.println(matrix.size()); System.out.println(matrix.get(4)); System.out.println(matrix.get(new int[]{0,0})); System.out.println(matrix.rowVector(1)); System.out.println(matrix.colVector(2)); Ranks<Integer> matrix1 = new Matrix<>(new Integer[][]{{7,8,9},{10,20,30}}); Addr<Integer> addr = new Addr<Integer>() { @Override public Integer zero() { return 0; } @Override public Integer one() { return 1; } @Override public Integer add(Integer lhs, Integer rhs) { return lhs + rhs; } @Override public Integer sub(Integer lhs, Integer rhs) { return lhs - rhs; } @Override public Integer mul(Integer k, Integer hs) { return k * hs; } @Override public Integer square(Integer hs) { return hs * hs; } @Override public double prescription(Integer hs) { return Math.sqrt(hs); } @Override public double div(Integer hs, double nhs) { return hs / nhs; } @Override public double acos(double hs) { return Math.acos(hs); } @Override public double toDegree(double hs) { return Math.toDegrees(hs); } }; System.out.println(matrix.trans()); Ranks<Integer> copy = ((Matrix)matrix).clone(); Attribute<Integer> vector = new Vector<>(new Integer[]{9,33,97}); Ranks<Integer> matrix2 = new Matrix<>(new Integer[][]{{7,9},{11,33},{42,97}}); //[(1,2,3),(3,4,5)]*(23,65,79) System.out.println(matrix.dot(vector,addr)); // [(1,2,3),(3,4,5)]*[(7,9),(11,33),(42,97)] System.out.println(matrix.dot(matrix2,addr)); Ranks<Integer> copy1 = ((Matrix)matrix).clone(); System.out.println(matrix.add(matrix1,addr)); System.out.println(matrix1.sub(copy,addr)); System.out.println(copy.div(2,addr)); System.out.println(copy.pos(addr,1)); System.out.println(copy.neg(addr,-1)); System.out.println(copy1.mul(3,addr)); System.out.println(Matrix.zero(2,3,addr)); Ranks<Integer> I = Matrix.identity(4, addr); System.out.println(I); Matrix<Integer> idenTest = new Matrix<>(new Integer[][]{{1,2,3,4},{5,6,7,8},{9,10,11,12},{13,14,15,16}}); System.out.println(idenTest.dot(I,addr)); System.out.println(I.dot(idenTest,addr)); } }
运行结果
Matrix{values=[[1, 2, 3],[3, 4, 5]]}
[2, 3]
2
3
6
4
1
Vector{values=[3, 4, 5]}
Vector{values=[3, 5]}
Matrix{values=[[1, 3],[2, 4],[3, 5]]}
Vector{values=[366, 644]}
Matrix{values=[[155, 366],[275, 644]]}
Matrix{values=[[8, 10, 12],[13, 24, 35]]}
Matrix{values=[[6, 6, 6],[7, 16, 25]]}
Matrix{values=[[0.5, 1.0, 1.5],[1.5, 2.0, 2.5]]}
Matrix{values=[[1, 2, 3],[3, 4, 5]]}
Matrix{values=[[-1, -2, -3],[-3, -4, -5]]}
Matrix{values=[[3, 6, 9],[9, 12, 15]]}
Matrix{values=[[0, 0, 0],[0, 0, 0]]}
Matrix{values=[[1, 0, 0, 0],[0, 1, 0, 0],[0, 0, 1, 0],[0, 0, 0, 1]]}
Matrix{values=[[1, 2, 3, 4],[5, 6, 7, 8],[9, 10, 11, 12],[13, 14, 15, 16]]}
Matrix{values=[[1, 2, 3, 4],[5, 6, 7, 8],[9, 10, 11, 12],[13, 14, 15, 16]]}
Python代码
from .Vector import Vector class Matrix: def __init__(self, list2d): self._values = [row[:] for row in list2d] @classmethod def zero(cls, r, c): # 返回一个r行c列的零矩阵 return cls([[0] * c for _ in range(r)]) @classmethod def identity(cls, n): # 返回一个n行n列的单位矩阵 m = [[0] * n for _ in range(n)] for i in range(n): m[i][i] = 1 return cls(m) def trans(self): # 返回矩阵的转置矩阵 return Matrix([[e for e in self.col_vector(i)] for i in range(self.col_num())]) def __add__(self, another): # 返回两个矩阵的加法结果 assert self.shape() == another.shape(), \ "加法错误，两个矩阵的形状必须相同" return Matrix([[a + b for a, b in zip(self.row_vector(i), another.row_vector(i))] for i in range(self.row_num())]) def __sub__(self, another): # 返回两个矩阵的减法结果 assert self.shape() == another.shape(), \ "减法错误，两个矩阵的形状必须相同" return Matrix([[a - b for a, b in zip(self.row_vector(i), another.row_vector(i))] for i in range(self.row_num())]) def dot(self, another): if isinstance(another, Vector): # 矩阵和向量的乘法 assert self.col_num() == len(another), \ "矩阵列数需要与向量的长度相等" return Vector([self.row_vector(i).dot(another) for i in range(self.row_num())]) if isinstance(another, Matrix): # 矩阵和矩阵的乘法 assert self.col_num() == another.row_num(), \ "当前矩阵列数必须等于点乘矩阵行数" return Matrix([[self.row_vector(i).dot(another.col_vector(j)) for j in range(another.col_num())] for i in range(self.row_num())]) def __mul__(self, k): # 返回矩阵的数量乘结果: self * k return Matrix([[e * k for e in self.row_vector(i)] for i in range(self.row_num())]) def __rmul__(self, k): # 返回矩阵的数量乘结果: k * self return self * k def __truediv__(self, k): # 返回数量除法的结果矩阵: self / k return (1 / k) * self def __pos__(self): # 返回矩阵取正的结果 return 1 * self def __neg__(self): # 返回矩阵取负的结果 return -1 * self def row_vector(self, index): # 返回矩阵的第index个行向量 return Vector(self._values[index]) def col_vector(self, index): # 返回矩阵的第index个列向量 return Vector([row[index] for row in self._values]) def __getitem__(self, pos): # 返回矩阵pos位置的元素 r, c = pos return self._values[r][c] def size(self): # 返回矩阵中元素的个数 r, c = self.shape() return r * c def row_num(self): # 返回矩阵的行数 return self.shape()[0] __len__ = row_num def col_num(self): # 返回矩阵的列数 return self.shape()[1] def shape(self): # 返回矩阵的形状：(行数，列数) return len(self._values), len(self._values[0]) def __repr__(self): return "Matrix({})".format(self._values) __str__ = __repr__
from playLA.Matrix import Matrix from playLA.Vector import Vector if __name__ == "__main__": matrix = Matrix([[1, 2, 3], [3, 4, 5]]) print(matrix) print("matrix.shape = {}".format(matrix.shape())) print("matrix.size = {}".format(matrix.size())) print("len(matrix) = {}".format(len(matrix))) print("matrix[0][0] = {}".format(matrix[0, 0])) print(matrix.row_vector(1)) print(matrix.col_vector(0)) matrix2 = Matrix([[7, 8, 9], [10, 20, 30]]) print("add: {}".format(matrix + matrix2)) print("subtract: {}".format(matrix2 - matrix)) print("scalar-mul: {}".format(matrix * 3)) print("scalar-mul: {}".format(3 * matrix)) print("zero_2_3: {}".format(Matrix.zero(2, 3))) vec = Vector([23, 65, 79]) print("matrix.dot(vec) = {}".format(matrix.dot(vec))) matrix3 = Matrix([[7, 9], [11, 33], [42, 97]]) print("matrix.dot(matrix3) = {}".format(matrix.dot(matrix3))) print("matrix.trans = {}".format(matrix.trans())) I = Matrix.identity(4) print(I) idenTest = Matrix([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]]) print(idenTest.dot(I)) print(I.dot(idenTest))
运行结果
Matrix([[1, 2, 3], [3, 4, 5]])
matrix.shape = (2, 3)
matrix.size = 6
len(matrix) = 2
matrix[0][0] = 1
(3, 4, 5)
(1, 3)
add: Matrix([[8, 10, 12], [13, 24, 35]])
subtract: Matrix([[6, 6, 6], [7, 16, 25]])
scalar-mul: Matrix([[3, 6, 9], [9, 12, 15]])
scalar-mul: Matrix([[3, 6, 9], [9, 12, 15]])
zero_2_3: Matrix([[0, 0, 0], [0, 0, 0]])
matrix.dot(vec) = (390, 724)
matrix.dot(matrix3) = Matrix([[155, 366], [275, 644]])
matrix.trans = Matrix([[1, 3], [2, 4], [3, 5]])
Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])
Matrix([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]])
Matrix([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]])
用矩阵表示空间
我们之前在矩阵与向量的点乘中，知道
矩阵与向量相乘，就是该矩阵的行向量与该向量的点乘，得到一个新的向量。例如
这其实是从行视角来看待的。但其实我们也可以通过列视角来看待矩阵与向量的点乘
它们二者其实是统一的，是同一个结果，只不过表达出来的形式是不同的。用图示来表示，即为
之前我们知道单位矩阵乘以一个矩阵依然为该矩阵本身，单位矩阵乘以向量依然如此，如
我们在之前的单位向量中知道，单位向量就是模为1的向量，在二维空间中，有两个标准单位向量
则上面的单位矩阵与向量相乘又可以表述为
，这就是(x,y)这个点坐标的定义！在(1,0)轴上的分量为x；在(0,1)轴上的分量为y
这就是通过列向量视角
最终得到的结果。
由于单位向量只表示方向，则我们可以把该单位矩阵与向量点乘的式子看成是定义出的空间坐标系的延展，其中
代表一个方向，
代表一个方向共建了一个平面空间。由此可见，矩阵定义了两个坐标轴。现在我们以任意一个矩阵
来看待这个问题，同理，它也有两个列向量
，它同样可以交织构成一个平面空间，在平面直角坐标系中
是这样的。
现在我们用该矩阵点乘一个向量
，则意味着矩阵中的两个列向量都沿着各自的方向放大两倍，最终得到
那么我们如何来理解这两个列向量交织的新的空间呢
我们将直角坐标系去掉，就变成了
而(2,2)这个点就是这个这个新空间坐标系中的红点
再把新空间的坐标转换成直角坐标系中的坐标。而之前的图形变换亦是如此，我们将变换矩阵想象成一个新的空间坐标系，就是将原图放入新的空间坐标系中的样子，只不过原图在原坐标系中的坐标原样不动的放入新空间坐标系中同样的坐标
比如图形沿着y轴翻转，它的变换矩阵为
，该变换矩阵所组成的新的空间坐标系如下
则我们把原图中的坐标原封不动的放入新的坐标系就变成了这样
而错切的变换矩阵为
，则它所交织的新的空间坐标系如下
同样，我们将原图坐标原封不动的放入新的坐标系中，就是这样
线性系统
什么是线性？未知数只能是一次方项，例如如下就是一个线性的方程
非线性方程：
因为这些方程在坐标系中的图像都不是直线，所以它们都不是线性的
而这样的方程得到的是一根直线，所以称为线性方程。
在二维的空间中，表现的是一根直线，而在三维的空间中则是一个平面
消元法
一个方程的左右两边同时乘以一个常数，一个方程加(减)另一个方程，交换两个方程的位置
高斯消元法
假设有这样一个方程组，根据之前的内容，我们可以知道，它可以等价于系数矩阵与未知数组成的向量相乘等于右边的常数向量
现在我们可以扔掉所有的未知数，而将常数向量也放入矩阵中
这样就形成了一个新的矩阵，我们称之为增广矩阵，我们所有的解方程的消元法的操作都可以在增广矩阵中完成。
首先我们可以用矩阵的第二个行向量减去第一个行向量乘以3，得到如下结果
再使用矩阵的第三个行向量减去第一个行向量乘以2，得到如下结果
再把矩阵的第三个行向量与第二个行向量相加，消去第三个行向的-1，得到如下结果
将矩阵的第三个行向量除以-15，得到结果如下
通过以上的过程，我们可以看出我们是在不断的消除第一个列向量中的元，我们称该元为主元
同样，第二个行向量的第二个元素，是第二行的主元，第三个行向量的第三个元素，是第三行的主元
所以我们每次消元的过程要首先确定每一行的主元是谁，就目前来看，其实第i行的主元就是第i个元素。而这种消元的方法就叫做高斯消元法。
但如果对于一个线性方程组的某一行没有主元该如何呢？
这里该方程组的增广矩阵为
由于第一个行向量没有主元，所以我们需要找到所有行中第一个元素最大的那个行向量交换位置，主要是为了尽量小的避免误差，就得到如下的形式，再去进行高斯消元法消元
再回来之前的例子，我们可以看到消元的目的是为了把所有的主元都变成1，而主元下面所有的元素都为0
现在这个增广矩阵实际上就相当于这样一个线性方程组
现在我们只要将所有行的除主元外的其他元素都变成0，就得到了该方程组的解。我们将矩阵的第二个行向量加上第三个行向量乘以10，就得到
此时我们将第一个行向量减去第三个行向量乘以4，得到
之后再将第一个行向量减去第二个行向量乘以2，得到
最终相当于
我们整个过程可以分为两部分来进行，第一部分为前向过程(从上到下)
- 选择最上的主元，化为1
- 主元下面的所有行减去主元所在行的某个倍数，使得主元下面所有元素都为0
第二部分为后向过程(从下到上)
- 选择最下的主元
- 主元上面的所有行减去主元所在行的某个倍数，使得主元上面的所有元素都为0
这两部分统称为高斯-约旦消元法，高斯负责了第一部分，约旦负责了第二部分
Addr接口增加相同类型的除法和比较大小的方法
public interface Addr<T extends Number> { T zero(); T one(); T add(T lhs,T rhs); T sub(T lhs,T rhs); T mul(T k,T hs); T square(T hs); double prescription(T hs); double div(T hs,double nhs); T div(T hs,T nhs); double acos(double hs); double toDegree(double hs); int compair(T lhs,T rhs); }
Ranks接口增加set方法
public interface Ranks<T extends Number> { /** * 矩阵的行列数 * @return */ int[] shape(); /** * 矩阵的行数 * @return */ int rowNum(); /** * 矩阵的列数 * @return */ int colNum(); /** * 矩阵的元素个数 * @return */ int size(); /** * 获取矩阵中的一个元素 * @param pos * @return */ T get(int pos); /** * 获取矩阵中的一个元素 * @param pos * @return */ T get(int[] pos); /** * 设置矩阵中某个位置的元素 * @param pos * @param value */ void set(int[] pos,T value); /** * 获取一个行向量 * @param index * @return */ Attribute<T> rowVector(int index); /** * 获取一个列向量 * @param index * @return */ Attribute<T> colVector(int index); /** * 矩阵相加 * @param another * @return */ Ranks<T> add(Ranks<T> another,Addr<T> addr); /** * 矩阵相减 * @param another * @param addr * @return */ Ranks<T> sub(Ranks<T> another,Addr<T> addr); /** * 矩阵乘以一个数 * @param k * @param addr * @return */ Ranks<T> mul(T k,Addr<T> addr); /** * 矩阵除以一个数 * @param k * @param addr * @return */ Ranks<Double> div(double k,Addr<T> addr); /** * 矩阵取正的结果 * @return */ Ranks<T> pos(Addr<T> addr,T one); /** * 矩阵取负的结果 * @return */ Ranks<T> neg(Addr<T> addr,T negOne); /** * 矩阵与向量的点乘 * @param another * @param addr * @return */ Attribute<T> dot(Attribute<T> another,Addr<T> addr); /** * 矩阵与矩阵的点乘 * @param another * @param addr * @return */ Ranks<T> dot(Ranks<T> another,Addr<T> addr); /** * 转置矩阵 * @return */ Ranks<T> trans(); }
矩阵实现类
public class Matrix<T extends Number> implements Ranks<T>,Cloneable { private T[][] values; public Matrix(T[][] values) { this.values = values; } @Override public String toString() { StringBuilder builder = new StringBuilder(); for (int i = 0; i < values.length - 1; i++) { builder.append(Arrays.toString(values[i]) + ","); } builder.append(Arrays.toString(values[values.length - 1])); return "Matrix{" + "values=" + "[" + builder.toString() + "]" + '}'; } /** * 创建一个零矩阵 * @param row * @param col * @param addr * @param <T> * @return */ @SuppressWarnings("unchecked") public static <T extends Number> Ranks<T> zero(int row,int col,Addr<T> addr) { T[][] values = (T[][])new Number[row][col]; for (int i = 0; i < row; i++) { Arrays.fill(values[i],addr.zero()); } return new Matrix<>(values); } /** * 创建一个单位矩阵 * @param n * @param <T> * @return */ @SuppressWarnings("unchecked") public static <T extends Number> Ranks<T> identity(int n,Addr<T> addr) { T[][] values = (T[][])new Number[n][n]; for (int i = 0; i < n; i++) { Arrays.fill(values[i],addr.zero()); } for (int i = 0; i < n; i++) { values[i][i] = addr.one(); } return new Matrix<>(values); } @Override public int[] shape() { int[] res = new int[2]; res[0] = values.length; res[1] = values[0].length; return res; } @Override public int rowNum() { return shape()[0]; } @Override public int colNum() { return shape()[1]; } @Override public int size() { int row = rowNum(); int col = colNum(); return row * col; } @Override public T get(int pos) { if (pos >= size() || pos < 0) { throw new IllegalArgumentException("位置超出索引范围"); } int row = pos / colNum(); int col = pos % colNum(); return values[row][col]; } @Override public T get(int[] pos) { if (pos.length != 2) { throw new IllegalArgumentException("位置长度超出范围"); } if (pos[0] >= rowNum() || pos[0] < 0 || pos[1] >= colNum() || pos[1] < 0) { throw new IllegalArgumentException("位置超出索引范围"); } return values[pos[0]][pos[1]]; } @Override public void set(int[] pos, T value) { if (pos.length != 2) { throw new IllegalArgumentException("位置长度超出范围"); } if (pos[0] >= rowNum() || pos[0] < 0 || pos[1] >= colNum() || pos[1] < 0) { throw new IllegalArgumentException("位置超出索引范围"); } values[pos[0]][pos[1]] = value; } @Override public Attribute<T> rowVector(int index) { if (index >= rowNum() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } return new Vector<>(values[index]); } @Override @SuppressWarnings("unchecked") public Attribute<T> colVector(int index) { if (index >= colNum() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } T[] values = (T[])new Number[rowNum()]; for (int i = 0; i < rowNum(); i++) { values[i] = this.values[i][index]; } return new Vector<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> add(Ranks<T> another,Addr<T> addr) { if (!Arrays.equals(this.shape(),another.shape())) { throw new IllegalArgumentException("加法错误，两个矩阵的形状必须相同"); } T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.add(this.get(new int[]{i,j}),another.get(new int[]{i,j})); } } this.values = values; return this; } @Override @SuppressWarnings("unchecked") public Ranks<T> sub(Ranks<T> another, Addr<T> addr) { if (!Arrays.equals(this.shape(),another.shape())) { throw new IllegalArgumentException("减法错误，两个矩阵的形状必须相同"); } T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.sub(this.get(new int[]{i,j}),another.get(new int[]{i,j})); } } this.values = values; return this; } @Override @SuppressWarnings("unchecked") public Ranks<T> mul(T k, Addr<T> addr) { T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.mul(k,this.get(new int[]{i,j})); } } this.values = values; return this; } @Override public Ranks<Double> div(double k, Addr<T> addr) { Double[][] values = new Double[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.div(this.get(new int[]{i,j}),k); } } return new Matrix<>(values); } @Override public Ranks<T> pos(Addr<T> addr,T one) { return mul(one,addr); } @Override public Ranks<T> neg(Addr<T> addr,T negOne) { return mul(negOne,addr); } @Override @SuppressWarnings("unchecked") public Attribute<T> dot(Attribute<T> another, Addr<T> addr) { if (this.colNum() != another.len()) { throw new IllegalArgumentException("矩阵列数需要与向量的长度相等"); } T[] values = (T[])new Number[this.rowNum()]; for (int i = 0; i < this.rowNum(); i++) { values[i] = this.rowVector(i).dot(another,addr); } return new Vector<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> dot(Ranks<T> another, Addr<T> addr) { if (this.colNum() != another.rowNum()) { throw new IllegalArgumentException("当前矩阵列数必须等于点乘矩阵行数"); } T[][] values = (T[][]) new Number[this.rowNum()][another.colNum()]; for (int i = 0; i < this.rowNum(); i++) { for (int j = 0; j < another.colNum(); j++) { values[i][j] = this.rowVector(i).dot(another.colVector(j),addr); } } return new Matrix<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> trans() { T[][] values = (T[][])new Number[this.colNum()][this.rowNum()]; for (int i = 0; i < this.colNum(); i++) { for (int j = 0; j < this.rowNum(); j++) { values[i][j] = this.values[j][i]; } } return new Matrix<>(values); } @Override @SuppressWarnings("unchecked") public Matrix clone() throws CloneNotSupportedException { T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = this.values[i][j]; } } return new Matrix(values); } }
新增一个接口
public interface FuncGroup<T extends Number> { /** * 高斯-约旦消元法 * @param addr */ void gaussJordanElimination(Addr<T> addr); /** * 打印增广矩阵 */ void fancyPrint(); }
线性系统实现类
/** * 线性系统 * @param <T> */ public class LinearSystem<T extends Number> implements FuncGroup<T>{ //系数矩阵的行数 private int matrixRow; //系数矩阵的列数 private int matrixCol; //增广矩阵 @Getter private Ranks<T> Ab; @SuppressWarnings("unchecked") public LinearSystem(Ranks<T> A,Attribute<T> b) { if (A.rowNum() != b.len()) { throw new IllegalArgumentException("矩阵的行数必须等于向量的长度"); } matrixRow = A.rowNum(); matrixCol = A.colNum(); if (matrixRow != matrixCol) { throw new IllegalArgumentException("要保证线性方程组具有唯一解"); } T[][] values = (T[][])new Number[A.rowNum()][A.colNum() + 1]; for (int i = 0; i < A.rowNum(); i++) { for (int j = 0; j <= A.colNum(); j++) { if (j < A.colNum()) { values[i][j] = A.get(new int[]{i, j}); }else { values[i][j] = b.get(i); } } } Ab = new Matrix<>(values); } @Override public void gaussJordanElimination(Addr<T> addr) { forward(addr); backword(addr); } @Override public void fancyPrint() { StringBuilder builder = new StringBuilder(); for (int i = 0; i < matrixRow; i++) { for (int j = 0; j < matrixCol + 1; j++) { if (j < matrixCol) { builder.append(" "); builder.append(Ab.get(new int[]{i, j})); } if (j == matrixCol) { builder.append(" | "); builder.append(Ab.get(new int[]{i,j})); builder.append("\n"); } } } System.out.println(builder.toString()); } /** * 前向 * @param addr */ @SuppressWarnings("unchecked") private void forward(Addr<T> addr) { int n = matrixRow; for (int i = 0; i < n; i++) { //Ab[i][i]为主元 int maxRow = maxRow(i, n, addr); T[] temp = (T[]) new Number[matrixCol + 1]; //每一行无论是否有主元，都与主元最大值所在行进行交换 for (int j = 0; j < matrixCol + 1; j++) { temp[j] = Ab.get(new int[]{i, j}); Ab.set(new int[]{i, j}, Ab.get(new int[]{maxRow, j})); Ab.set(new int[]{maxRow, j}, temp[j]); } T u = Ab.get(new int[]{i, i}); //将主元归为一 for (int j = 0; j < matrixCol + 1; j++) { Ab.set(new int[]{i, j}, addr.div(Ab.get(new int[]{i, j}), u)); } //将主元所在行以下的行的所有非主元数据全部变成0 for (int k = i + 1; k < n; k++) { T u1 = Ab.get(new int[]{k, i}); for (int l = 0; l < matrixCol + 1; l++) { Ab.set(new int[]{k, l}, addr.sub(Ab.get(new int[]{k, l}), addr.mul(u1, Ab.get(new int[]{i, l})))); } } } } /** * 后向 */ private void backword(Addr<T> addr) { int n = matrixRow; for (int i = n - 1; i >= 0; i--) { //Ab[i][i]为主元 for (int j = i - 1; j >= 0; j--) { T u = Ab.get(new int[]{j, i}); for (int k = 0; k < matrixCol + 1; k++) { Ab.set(new int[]{j,k},addr.sub(Ab.get(new int[]{j,k}), addr.mul(u,Ab.get(new int[]{i,k})))); } } } } /** * 获取主元最大值所在的行 * @param index * @param n * @param addr * @return */ private int maxRow(int index,int n,Addr<T> addr) { T best = Ab.get(new int[]{index,index}); int ret = index; for (int i = index + 1; i < n; i++) { if (addr.compair(Ab.get(new int[]{i,index}),best) == 1) { best = Ab.get(new int[]{i,index}); ret = i; } } return ret; } public static void main(String[] args) { Ranks<Double> A = new Matrix<>(new Double[][]{{1.0,2.0,4.0},{3.0,7.0,2.0},{2.0,3.0,3.0}}); Attribute<Double> b = new Vector<>(new Double[]{7.0,-11.0,1.0}); FuncGroup<Double> ls = new LinearSystem<>(A,b); Addr<Double> addr = new Addr<Double>() { @Override public Double zero() { return 0.0; } @Override public Double one() { return 1.0; } @Override public Double add(Double lhs, Double rhs) { return lhs + rhs; } @Override public Double sub(Double lhs, Double rhs) { return lhs - rhs; } @Override public Double mul(Double k, Double hs) { return k * hs; } @Override public Double square(Double hs) { return hs * hs; } @Override public double prescription(Double hs) { return Math.sqrt(hs); } @Override public double div(Double hs, double nhs) { return hs / nhs; } @Override public Double div(Double hs, Double nhs) { return hs / nhs; } @Override public double acos(double hs) { return Math.acos(hs); } @Override public double toDegree(double hs) { return Math.toDegrees(hs); } @Override public int compair(Double lhs, Double rhs) { if (lhs - rhs > 0) { return 1; }else if (lhs - rhs < 0) { return -1; }else { return 0; } } }; ls.gaussJordanElimination(addr); ls.fancyPrint(); } }
运行结果
1.0 0.0 0.0 | -0.9999999999999991
-0.0 1.0 0.0 | -2.0
-0.0 -0.0 1.0 | 2.9999999999999996
Python代码
向量类增加部分方法
import math class Vector: def __init__(self, lst): self._values = lst @classmethod def zero(cls,dim): # 返回一个dim维的零向量 return cls([0] * dim) def project(self, another): # 其他向量在当前向量的投影 return self.normalize() * (self.dot(another) / self.norm()) def degree(self, another): # 返回向量的夹角 return math.degrees(math.acos(self.dot(another) / (self.norm() * another.norm()))) def norm(self): # 返回向量的模 return math.sqrt(sum(e**2 for e in self)) def normalize(self): # 返回向量的单位向量 return Vector([e / self.norm() for e in self]) def underlying_list(self): # 返回向量的底层列表 return self._values[:] def __add__(self, another): # 向量加法，返回结果向量 assert len(self) == len(another),\ "加法错误，两个向量的维度必须相等" return Vector([a + b for a, b in zip(self, another)]) def __sub__(self, another): # 向量减法，返回结果向量 assert len(self) == len(another), \ "减法错误，两个向量的维度必须相等" return Vector([a - b for a, b in zip(self, another)]) def __pos__(self): # 返回向量取正的结果向量 return 1 * self def __neg__(self): # 返回向量取负的结果向量 return -1 * self def dot(self, another): # 向量点乘，返回结果标量 assert len(self) == len(another), \ "点乘错误，两个向量的维度必须相等" return sum(a * b for a, b in zip(self, another)) def __mul__(self, k): # 返回数量乘法的结果向量:self * k return Vector([k * e for e in self]) def __rmul__(self, k): # 返回数量乘法的结果向量:k * self return self * k def __truediv__(self, k): # 返回数量除法的结果向量:self / k return (1 / k) * self def __iter__(self): # 返回向量的迭代器 return self._values.__iter__() def __getitem__(self, index): # 取向量的第index个元素 return self._values[index] def __len__(self): # 返回向量的长度 return len(self._values) def __repr__(self): return "Vector({})".format(self._values) def __str__(self): return "({})".format(", ".join(str(e) for e in self._values))
线性系统类
from .Vector import Vector class LinearSystem: def __init__(self, A, b): assert A.row_num() == len(b), "矩阵的行数必须等于向量的长度" self._matrix_row = A.row_num() self._matrix_col = A.col_num() assert self._matrix_row == self._matrix_col, "要保证线性方程组具有唯一解" self.Ab = [Vector(A.row_vector(i).underlying_list() + [b[i]]) for i in range(self._matrix_row)] def _max_row(self, index, n): # 主元最大值所在的行 best, ret = self.Ab[index][index], index for i in range(index + 1, n): if self.Ab[i][index] > best: best, ret = self.Ab[i][index], i return ret def _forward(self): # 前向 n = self._matrix_row for i in range(n): # Ab[i][i]为主元 max_row = self._max_row(i, n) self.Ab[i], self.Ab[max_row] = self.Ab[max_row], self.Ab[i] # 将主元归为一 self.Ab[i] = self.Ab[i] / self.Ab[i][i] for j in range(i + 1, n): self.Ab[j] = self.Ab[j] - self.Ab[j][i] * self.Ab[i] def _backward(self): # 后向 n = self._matrix_row for i in range(n - 1, -1, -1): # Ab[i][i]为主元 for j in range(i - 1, -1, -1): self.Ab[j] = self.Ab[j] - self.Ab[j][i] * self.Ab[i] def gauss_jordan_elimination(self): # 高斯-约旦消元法 self._forward() self._backward() def fancy_print(self): # 打印增广矩阵 for i in range(self._matrix_row): print(" ".join(str(self.Ab[i][j]) for j in range(self._matrix_col)), end=" ") print("|", self.Ab[i][-1])
from playLA.Matrix import Matrix from playLA.Vector import Vector from playLA.LinearSystem import LinearSystem if __name__ == "__main__": A = Matrix([[1, 2, 4], [3, 7, 2], [2, 3, 3]]) b = Vector([7, -11, 1]) ls = LinearSystem(A, b) ls.gauss_jordan_elimination() ls.fancy_print()
运行结果
1.0 0.0 0.0 | -1.0
-0.0 1.0 0.0 | -2.0
-0.0 -0.0 1.0 | 3.0
行最简形式和线性方程组解的结构
之前我们在代码中可以看到这样一段
if (matrixRow != matrixCol) { throw new IllegalArgumentException("要保证线性方程组具有唯一解"); }
但一个线性方程组未必具有唯一解，它可能有无数个解，也可能没有解。
我们来看这样一组线性方程组
则它的增广矩阵如下
我们用第二个行向量加上第一个行向量，第三个行向量减去第一个行向量乘以2得
再用第三个行向量加上第二个行向量，得
现在第三行表示
，这显然是不可能的，此时这个线性方程组无解
我们再来看一组线性方程组
则它的增广矩阵如下
我们将第二个行向量加上第一个行向量，第三个行向量减去第一个行向量乘以3，得
我们将第二行进行归一操作，得
现在我们将第三个行向量加上第二个行向量，得
再将第一个行向量减去第二个行向量乘以2，得
现在这个增广矩阵表示的是
经过转化，最后变为
这就是该线性系统的解，z任意取值，都能得到一组x,y,z满足方程组，则方程组有无数组解。
无论一个线性方程组有唯一解，无解还是有无数组解，它们都满足一种阶梯型，称为阶梯型矩阵
如上面的三种情况，我们都可以绘制成阶梯形状
根据上图，我们能很清楚的看到，在阶梯的下层的元素都为0。而阶梯的上层与阶梯相邻的都是它的主元。非零行的第一个元素(主元)为1.主元所在列的其他元素均为0.而经过了高斯-约旦消元法之后得到的这个矩阵，我们称为——行最简形式(reduced row echelon form(RREF))
对于第一行来说，主元未必是第一个元素，而在下面的第i行来说，主元的位置也未必在第i列上
而以下的例子中则不是行最简形式
现在我们可以来总结一下方程组的解的各种情况
而满足上面三种情况，需要达到下面的条件
- 方程组有唯一解
系数矩阵非零行=未知数个数
- 方程组无解
系数矩阵非零行<行最简形式非零行(增广矩阵非零行)
- 方程组有无数解
系数矩阵非零行<未知数个数
在下面的增广矩阵中
我们可以判断出，该方程组是有无数解的，因为系数矩阵非零行3<未知数个数4
而
是有唯一解的，系数矩阵非零行3等于未知数个数3
更一般的线性系统求解
面对上面的不同解的情况，我们就需要调整方法。而高斯-约旦消元法的主要思想不变，但需要调整细节。
这是一个有四个方程，五个未知数的方程组，转化为增广矩阵如下
使用高斯-约旦消元法，我们将第二个行向量与第一个行向量相加，第三个行向量减去第一个行向量，第四个行向量加上第一个行向量乘以2，得到
这里第一行第一个1即为首元。这里我们可以看到，第一行下面，不仅是第一列变成0了，第二列也变成0了。按照之前的情况，如果第二行第二列的元素为0，则我们会向下找一个不为0的最大的元素进行行交换，来获取第二行的主元。但是很可惜，第二行之下第二列之下，所有元素都为0。此时就需要继续相后寻找第二行的主元，那就是第三列的2作为第二行的主元，再将该主元归一，得到
再将该主元下面的元素变成0.第三个行向量减去第二个行向量乘以2，第四个行向量加上第二个行向量，得到
跟之前的情况一样，我们在第三行第四列也找不到主元，于是我们继续向后寻找主元，便找到第三行第五列的3作为主元。再将该主元归一。
将第三行主元下面的元素变成0.第四个行向量减去第三个行向量乘以2，得到
我们可以看到，第四行已经全部变成了0，则第四行没有主元。现在我们要将最后一个主元上面的元素变成0，第一个行向量减去第三个行向量乘以3，第二个行向量加上第三个行向量，得到
此时我们需要将第二行的主元上面的元素化为0.第一个行向量减去第二个行向量乘以2，得到
此时我们的高斯-约旦消元法到此就结束了，得到了行最简形式。由之前的系数矩阵非零行3<未知数个数5可知，该方程组有无数个解
跟唯一解的线性系统不同，它不是所有的列都包含了主元，我们把含有主元的列叫做主元列，如这里的第一列、第三列、第五列为主元列。而不含有主元的列，称为自由列，如这里的第二列、第四列为自由列。而把该行最简形式化成方程组，则为
其中x、z、u是主元列的未知数，而y、w是自由列的未知数，而自由列的未知数可以任意取值，得到方程组的无数个解。
Addr接口增加判断两个元素相等的方法
public interface Addr<T extends Number> { T zero(); T one(); T add(T lhs,T rhs); T sub(T lhs,T rhs); T mul(T k,T hs); T square(T hs); double prescription(T hs); double div(T hs,double nhs); T div(T hs,T nhs); double acos(double hs); double toDegree(double hs); int compair(T lhs,T rhs); boolean equals(T lhs,T rhs); }
FuncGroup接口修改高斯-约旦消元法方法，返回是否有解
public interface FuncGroup<T extends Number> { /** * 高斯-约旦消元法，并判断是否有解 * @param addr */ boolean gaussJordanElimination(Addr<T> addr); /** * 打印增广矩阵 */ void fancyPrint(); }
线性系统实现类
/** * 线性系统 * @param <T> */ public class LinearSystem<T extends Number> implements FuncGroup<T>{ //系数矩阵的行数 private int matrixRow; //系数矩阵的列数 private int matrixCol; //增广矩阵 @Getter private Ranks<T> Ab; //主元列 @Getter private List<Integer> pivots = new ArrayList<>(); @SuppressWarnings("unchecked") public LinearSystem(Ranks<T> A,Attribute<T> b) { if (A.rowNum() != b.len()) { throw new IllegalArgumentException("矩阵的行数必须等于向量的长度"); } matrixRow = A.rowNum(); matrixCol = A.colNum(); // if (matrixRow != matrixCol) { // throw new IllegalArgumentException("要保证线性方程组具有唯一解"); // } T[][] values = (T[][])new Number[A.rowNum()][A.colNum() + 1]; for (int i = 0; i < A.rowNum(); i++) { for (int j = 0; j <= A.colNum(); j++) { if (j < A.colNum()) { values[i][j] = A.get(new int[]{i, j}); }else { values[i][j] = b.get(i); } } } Ab = new Matrix<>(values); } @Override public boolean gaussJordanElimination(Addr<T> addr) { forward(addr); backword(addr); //无主元的行的非系数元素如果为0则无解 for (int i = pivots.size(); i < matrixRow; i++) { if (!addr.equals(Ab.get(new int[]{i,matrixCol}),addr.zero())) { return false; } } return true; } @Override public void fancyPrint() { StringBuilder builder = new StringBuilder(); for (int i = 0; i < matrixRow; i++) { for (int j = 0; j < matrixCol + 1; j++) { if (j < matrixCol) { builder.append(" "); builder.append(Ab.get(new int[]{i, j})); } if (j == matrixCol) { builder.append(" | "); builder.append(Ab.get(new int[]{i,j})); builder.append("\n"); } } } System.out.println(builder.toString()); } /** * 前向 * @param addr */ @SuppressWarnings("unchecked") private void forward(Addr<T> addr) { // int n = matrixRow; int i = 0; int k = 0; // for (int i = 0; i < n; i++) { while (i < matrixRow && k < matrixCol) { //看Ab[i][k]位置是否可以是主元 int maxRow = maxRow(i, k, matrixRow, addr); T[] temp = (T[]) new Number[matrixCol + 1]; //每一行无论是否有主元，都与主元最大值所在行进行交换 for (int j = 0; j < matrixCol + 1; j++) { temp[j] = Ab.get(new int[]{i, j}); Ab.set(new int[]{i, j}, Ab.get(new int[]{maxRow, j})); Ab.set(new int[]{maxRow, j}, temp[j]); } if (addr.equals(Ab.get(new int[]{i,k}),addr.zero())) { k++; }else { T u = Ab.get(new int[]{i, k}); //将主元归为一 for (int j = 0; j < matrixCol + 1; j++) { Ab.set(new int[]{i, j}, addr.div(Ab.get(new int[]{i, j}), u)); } //将主元所在行以下的行的所有非主元数据全部变成0 for (int m = i + 1; m < matrixRow; m++) { T u1 = Ab.get(new int[]{m, k}); for (int l = 0; l < matrixCol + 1; l++) { Ab.set(new int[]{m, l}, addr.sub(Ab.get(new int[]{m, l}), addr.mul(u1, Ab.get(new int[]{i, l})))); } } pivots.add(k); i++; } } } /** * 后向 */ private void backword(Addr<T> addr) { int n = pivots.size(); for (int i = n - 1; i >= 0; i--) { int k = pivots.get(i); //Ab[i][k]为主元 for (int j = i - 1; j >= 0; j--) { T u = Ab.get(new int[]{j, k}); for (int l = 0; l < matrixCol + 1; l++) { Ab.set(new int[]{j,l},addr.sub(Ab.get(new int[]{j,l}), addr.mul(u,Ab.get(new int[]{i,l})))); } } } } /** * 获取主元最大值所在的行 * @param addr * @return */ private int maxRow(int indexI,int indexJ,int n,Addr<T> addr) { T best = Ab.get(new int[]{indexI,indexJ}); int ret = indexI; for (int i = indexI + 1; i < n; i++) { if (addr.compair(Ab.get(new int[]{i,indexJ}),best) == 1) { best = Ab.get(new int[]{i,indexJ}); ret = i; } } return ret; } public static void main(String[] args) { Ranks<Double> A = new Matrix<>(new Double[][]{{1.0,2.0,4.0},{3.0,7.0,2.0},{2.0,3.0,3.0}}); Attribute<Double> b = new Vector<>(new Double[]{7.0,-11.0,1.0}); FuncGroup<Double> ls = new LinearSystem<>(A,b); Addr<Double> addr = new Addr<Double>() { @Override public Double zero() { return 0.0; } @Override public Double one() { return 1.0; } @Override public Double add(Double lhs, Double rhs) { return lhs + rhs; } @Override public Double sub(Double lhs, Double rhs) { return lhs - rhs; } @Override public Double mul(Double k, Double hs) { return k * hs; } @Override public Double square(Double hs) { return hs * hs; } @Override public double prescription(Double hs) { return Math.sqrt(hs); } @Override public double div(Double hs, double nhs) { return hs / nhs; } @Override public Double div(Double hs, Double nhs) { return hs / nhs; } @Override public double acos(double hs) { return Math.acos(hs); } @Override public double toDegree(double hs) { return Math.toDegrees(hs); } @Override public int compair(Double lhs, Double rhs) { if (lhs - rhs > 0) { return 1; }else if (lhs - rhs < 0) { return -1; }else { return 0; } } @Override public boolean equals(Double lhs, Double rhs) { double dis = 1e-8; if (Math.abs(lhs - rhs) < dis) { return true; } return false; } }; if (!ls.gaussJordanElimination(addr)) { System.out.println(" 无解"); } ls.fancyPrint(); System.out.println(); Ranks<Double> A1 = new Matrix<>(new Double[][]{{1.0,-1.0,2.0,0.0,3.0}, {-1.0,1.0,0.0,2.0,-5.0},{1.0,-1.0,4.0,2.0,4.0},{-2.0,2.0,-5.0,-1.0,-3.0}}); Attribute<Double> b1 = new Vector<>(new Double[]{1.0,5.0,13.0,-1.0}); FuncGroup<Double> ls1 = new LinearSystem<>(A1,b1); if (!ls1.gaussJordanElimination(addr)) { System.out.println(" 无解"); } ls1.fancyPrint(); System.out.println(); Ranks<Double> A2 = new Matrix<>(new Double[][]{{2.0,2.0},{2.0,1.0},{1.0,2.0}}); Attribute<Double> b2 = new Vector<>(new Double[]{3.0,2.5,7.0}); FuncGroup<Double> ls2 = new LinearSystem<>(A2,b2); if (!ls2.gaussJordanElimination(addr)) { System.out.println(" 无解"); } ls2.fancyPrint(); } }
运行结果
1.0 0.0 0.0 | -0.9999999999999991
-0.0 1.0 0.0 | -2.0
-0.0 -0.0 1.0 | 2.9999999999999996
1.0 -1.0 0.0 -2.0 0.0 | -15.0
0.0 0.0 1.0 1.0 0.0 | 5.0
0.0 0.0 0.0 0.0 1.0 | 2.0
0.0 0.0 0.0 0.0 0.0 | 0.0
无解
1.0 0.0 | -4.0
0.0 1.0 | 5.5
0.0 0.0 | 5.0
Python代码
EPSILON = 1e-8 def is_zero(x): return abs(x) < EPSILON def is_equal(a, b): return abs(a - b) < EPSILON
from .Vector import Vector from ._global import is_zero class LinearSystem: def __init__(self, A, b): assert A.row_num() == len(b), "矩阵的行数必须等于向量的长度" self._matrix_row = A.row_num() self._matrix_col = A.col_num() self.Ab = [Vector(A.row_vector(i).underlying_list() + [b[i]]) for i in range(self._matrix_row)] self.pivots = [] def _max_row(self, index_i, index_j, n): # 主元最大值所在的行 best, ret = self.Ab[index_i][index_j], index_i for i in range(index_i + 1, n): if self.Ab[i][index_j] > best: best, ret = self.Ab[i][index_j], i return ret def _forward(self): # 前向 i, k = 0, 0 while i < self._matrix_row and k < self._matrix_col: # 看Ab[i][k]位置是否可以为主元 max_row = self._max_row(i, k, self._matrix_row) self.Ab[i], self.Ab[max_row] = self.Ab[max_row], self.Ab[i] if is_zero(self.Ab[i][k]): k += 1 else: # 将主元归为一 self.Ab[i] = self.Ab[i] / self.Ab[i][k] for j in range(i + 1, self._matrix_row): self.Ab[j] = self.Ab[j] - self.Ab[j][k] * self.Ab[i] self.pivots.append(k) i += 1 def _backward(self): # 后向 n = len(self.pivots) for i in range(n - 1, -1, -1): k = self.pivots[i] # Ab[i][k]为主元 for j in range(i - 1, -1, -1): self.Ab[j] = self.Ab[j] - self.Ab[j][k] * self.Ab[i] def gauss_jordan_elimination(self): # 高斯-约旦消元法，如果有解，返回True;如果没有解，返回False self._forward() self._backward() for i in range(len(self.pivots), self._matrix_row): if not is_zero(self.Ab[i][-1]): return False return True def fancy_print(self): # 打印增广矩阵 for i in range(self._matrix_row): print(" ".join(str(self.Ab[i][j]) for j in range(self._matrix_col)), end=" ") print("|", self.Ab[i][-1])
from playLA.Matrix import Matrix from playLA.Vector import Vector from playLA.LinearSystem import LinearSystem if __name__ == "__main__": A = Matrix([[1, 2, 4], [3, 7, 2], [2, 3, 3]]) b = Vector([7, -11, 1]) ls = LinearSystem(A, b) ls.gauss_jordan_elimination() ls.fancy_print() print() A1 = Matrix([[1, -1, 2, 0, 3], [-1, 1, 0, 2, -5], [1, -1, 4, 2, 4], [-2, 2, -5, -1, -3]]) b1 = Vector([1, 5, 13, -1]) ls1 = LinearSystem(A1, b1) ls1.gauss_jordan_elimination() ls1.fancy_print() print() A2 = Matrix([[2, 2], [2, 1], [1, 2]]) b2 = Vector([3, 2.5, 7]) ls2 = LinearSystem(A2, b2) if not ls2.gauss_jordan_elimination(): print("无解") ls2.fancy_print()
运行结果
1.0 0.0 0.0 | -1.0
-0.0 1.0 0.0 | -2.0
-0.0 -0.0 1.0 | 3.0
1.0 -1.0 0.0 -2.0 0.0 | -15.0
0.0 0.0 1.0 1.0 0.0 | 5.0
0.0 0.0 0.0 0.0 1.0 | 2.0
0.0 0.0 0.0 0.0 0.0 | 0.0
无解
1.0 0.0 | -4.0
0.0 1.0 | 5.5
0.0 0.0 | 5.0
直观理解线性方程组解的结构
n个未知数有n个方程，才可能有唯一解
如果有一个方程为2x + 2y = 3,在直角坐标系中，它可以表示为一条直线，这条直线上的任意一个点都满足这个方程
这个时候我们再增加一个方程2x + y = 2，那么它也是一根直线
很显然，这两根直线相交的点就是方程组的唯一解
当然，在这种情况下它还有可能无解，如
如果有3个未知数，如x + y + 2z = 0，在三维空间中，它就是一个平面
如果有两个三元方程联立，则可能表示两个平面相交，它是一条直线
当方程组个数小于未知数个数，一定没有唯一解，而这两个平面相交的直线代表了无数个解。当然如果两个平面是平行的，则代表方程组无解
三个三元方程联立，它们代表三个平面，它们有可能相交于一个点
这个时候，方程组就有唯一解
它们也有可能相交于一条线
此时方程组就有无数解
如果三个平面平行
此时方程组就无解，但是三个三元方程联立没有解不仅仅是这一种情况
以上的情况，方程组都是无解的。
四个三元方程联立，情况也是类似的，代表了四个平面的各种情况
同理，三个二元方程联立，就代表了三条直线，它们有可能无解
也可能有唯一解
也可能有无数解(三个方程是同一条直线)
则方程个数和未知数个数的关系为
| 方程个数<未知数 | 方程个数=未知数 | 方程个数>未知数 |
| 无解 | 无解 | 无解 |
| 唯一解 | 唯一解 | |
| 无数解 | 无数解 | 无数解 |
对于行最简形式，系数矩阵的非零行不可能大于未知数个数
| 行最简形式系数矩阵非零行<未知数 | 行最简形式系数矩阵非零行=未知数 |
| 无数解 | 唯一解 |
以上是当方程组有解的情况
行最简形式系数矩阵非零行<行最简形式增广矩阵非零行，则无解
齐次线性方程组
齐次线性方程组的特性就是它的等号的右面都为0.满足这个条件的方程组就叫齐次线性方程组。
齐次线性方程和我们求解普通线性方程组的过程没有任何区别。首先划成增广矩阵
再第一行的主元归一，将主元下面的元素变成0
再将第二行的主元归一，将主元下面的元素变成0
最后化成行最简形式
对于齐次线性方程组来说，是一定有解的。因为对于齐次线性方程组来说，如果每一个未知数都取0，则肯定是它的一个解。对于齐次线性方程组来说，我们要看的是它是有唯一的0解还是它有无数个解。这个判断方式跟普通的线性方程组的方式是一样的。对于上面的这个行最简形式来说，它的系数矩阵的非零行的个数2小于它的未知数的个数3，所以它有无数个解。由于它的最后一列永远为0，所以它肯定不会出现普通线性方程组可能出现的那个矛盾——系数矩阵非零行<增广矩阵非零行，而导致的无解。
而我们在齐次线性方程组的计算中，也可以不对增广矩阵进行操作，只看它的系数矩阵就可以了。
当然我们也可以不这么操作，而是只当作普通线性方程组来看，使用增广矩阵来进行操作。
当然相对齐次线性方程组来说，就有非齐次线性方程组，那就是等号右边的值不都为0
矩阵的逆
在数字系统中，有
，而在矩阵中，如果有AB=BA=I(I是单位矩阵)，则称B是A的逆矩阵，记做
由于数字乘法满足乘法交换律，但在矩阵的乘法并不满足乘法交换律，如果A能找到逆矩阵的话，A称为可逆矩阵，或者叫做非奇异矩阵(non-singular)。有些矩阵是不可逆的，称为不可逆矩阵，或者奇异矩阵(singular)
比如有这样一个矩阵，则它的逆矩阵为
如果BA=I，则称B是A的左逆矩阵
如果AC=I, 则称C是A的右逆矩阵
如果一个矩阵A即存在左逆矩阵B,又存在右逆矩阵C,则B=C。对于矩阵A,存在矩阵B,满足BA=AB=I,矩阵A可逆
可逆矩阵一定是方阵,非方阵一定不可逆
现在我们知道了
，
就是矩阵的逆
求解矩阵的逆
我们还是以
来说明
由矩阵的乘法，由于A是一个2*2的矩阵，而I也是一个2*2的矩阵，则A的逆也一定是一个2*2的矩阵，则我们给A的逆矩阵设为
则满足矩阵乘法
则我们可以把该式子拆成
变换后为
则它们可以转变成这样一个线性系统
由于第一、第二个方程和第三、第四个方程的未知数并不想干，我们可以分别将第一、第二；第三、第四个方程组分别化成两个增广矩阵
将其通过高斯-约旦消元法化成行最简形式，可得
由于两个增广矩阵的系数矩阵相同，我们可以将其合并成一个新的增广矩阵
再通过高斯-约旦消元法化成行最简形式
则可以得出，对于任意一个有逆矩阵的矩阵A有如下的形式
现在我们用高斯-约旦消元法来求出这个逆矩阵
我们将第二个行向量减去第一个行向量乘以3，得到
再将第二行的主元归一
再将第一个行向量减去第二个行向量乘以2，得
至此，我们就化成了行最简形式，则我们可以看出原矩阵的逆矩阵就是
现在我们来看一个这个过程的一般性
首先我们可以看出，这个增广矩阵是不可能有无数解的，因为有无数解的前提是有一个化为行最简形式的时候有一个全零行，而该增广矩阵的右边是一个单位矩阵，无论你怎么变换，都不可能将单位矩阵的元素全部变成0.也就是说，对于一个方阵，要不然我们可以求出唯一一个矩阵的逆，要不然我们求不出矩阵的逆。如果一个矩阵存在逆矩阵，则它的逆矩阵是唯一的。但这个线性系统也可能是无解的，即没有逆矩阵，如果在变成行最简形式后，系数矩阵的某一行变成全零行，此时就是无解的，因为单位矩阵的每一行是必有值的，不可能变成全零行。
如果一个方阵A有右逆B,则B也是A的左逆，即B是A的逆。
给Ranks接口增加一个求矩阵的逆的方法
public interface Ranks<T extends Number> { /** * 矩阵的行列数 * @return */ int[] shape(); /** * 矩阵的行数 * @return */ int rowNum(); /** * 矩阵的列数 * @return */ int colNum(); /** * 矩阵的元素个数 * @return */ int size(); /** * 获取矩阵中的一个元素 * @param pos * @return */ T get(int pos); /** * 获取矩阵中的一个元素 * @param pos * @return */ T get(int[] pos); /** * 设置矩阵中某个位置的元素 * @param pos * @param value */ void set(int[] pos,T value); /** * 获取一个行向量 * @param index * @return */ Attribute<T> rowVector(int index); /** * 获取一个列向量 * @param index * @return */ Attribute<T> colVector(int index); /** * 矩阵相加 * @param another * @return */ Ranks<T> add(Ranks<T> another,Addr<T> addr); /** * 矩阵相减 * @param another * @param addr * @return */ Ranks<T> sub(Ranks<T> another,Addr<T> addr); /** * 矩阵乘以一个数 * @param k * @param addr * @return */ Ranks<T> mul(T k,Addr<T> addr); /** * 矩阵除以一个数 * @param k * @param addr * @return */ Ranks<Double> div(double k,Addr<T> addr); /** * 矩阵取正的结果 * @return */ Ranks<T> pos(Addr<T> addr,T one); /** * 矩阵取负的结果 * @return */ Ranks<T> neg(Addr<T> addr,T negOne); /** * 矩阵与向量的点乘 * @param another * @param addr * @return */ Attribute<T> dot(Attribute<T> another,Addr<T> addr); /** * 矩阵与矩阵的点乘 * @param another * @param addr * @return */ Ranks<T> dot(Ranks<T> another,Addr<T> addr); /** * 转置矩阵 * @return */ Ranks<T> trans(); /** * 矩阵的逆 * @param addr * @return */ Ranks<T> inv(Addr<T> addr); }
修改线性系统类，增加一个新的构造方法
/** * 线性系统 * @param <T> */ public class LinearSystem<T extends Number> implements FuncGroup<T>{ //系数矩阵的行数 private int matrixRow; //系数矩阵的列数 private int matrixCol; //增广矩阵 @Getter private Ranks<T> Ab; //主元列 @Getter private List<Integer> pivots = new ArrayList<>(); @SuppressWarnings("unchecked") public LinearSystem(Ranks<T> A,Attribute<T> b) { if (A.rowNum() != b.len()) { throw new IllegalArgumentException("系数矩阵的行数必须等于常数向量的长度"); } matrixRow = A.rowNum(); matrixCol = A.colNum(); T[][] values = (T[][])new Number[A.rowNum()][A.colNum() + 1]; for (int i = 0; i < A.rowNum(); i++) { for (int j = 0; j <= A.colNum(); j++) { if (j < A.colNum()) { values[i][j] = A.get(new int[]{i, j}); }else { values[i][j] = b.get(i); } } } Ab = new Matrix<>(values); } @SuppressWarnings("unchecked") public LinearSystem(Ranks<T> A,Ranks<T> b) { if (A.rowNum() != b.colNum()) { throw new IllegalArgumentException("系数矩阵的行数必须等于单位矩阵的列数"); } matrixRow = A.rowNum(); matrixCol = A.colNum(); T[][] values = (T[][])new Number[A.rowNum()][A.colNum() + b.colNum()]; for (int i = 0; i < A.rowNum(); i++) { for (int j = 0; j < A.colNum() + b.colNum(); j++) { if (j < A.colNum()) { values[i][j] = A.get(new int[]{i, j}); }else { values[i][j] = b.get(new int[]{i, j - A.colNum()}); } } } Ab = new Matrix<>(values); } @Override public boolean gaussJordanElimination(Addr<T> addr) { forward(addr); backword(addr); //无主元的行的非系数元素如果为0则无解 for (int i = pivots.size(); i < matrixRow; i++) { if (!addr.equals(Ab.get(new int[]{i,matrixCol}),addr.zero())) { return false; } } return true; } @Override public void fancyPrint() { StringBuilder builder = new StringBuilder(); if (Ab.colNum() - matrixCol > 1) { for (int i = 0; i < matrixRow; i++) { for (int j = 0; j < Ab.colNum(); j++) { if (j < matrixCol) { builder.append(" "); builder.append(Ab.get(new int[]{i, j})); } else if (j == matrixCol) { builder.append(" | "); builder.append(Ab.get(new int[]{i, j})); } else { builder.append(" "); builder.append(Ab.get(new int[]{i, j})); builder.append("\n"); } } } }else if (Ab.colNum() - matrixCol == 1) { for (int i = 0; i < matrixRow; i++) { for (int j = 0; j < Ab.colNum(); j++) { if (j < matrixCol) { builder.append(" "); builder.append(Ab.get(new int[]{i, j})); } else if (j == matrixCol) { builder.append(" | "); builder.append(Ab.get(new int[]{i, j})); builder.append("\n"); } } } } System.out.println(builder.toString()); } /** * 前向 * @param addr */ @SuppressWarnings("unchecked") private void forward(Addr<T> addr) { // int n = matrixRow; int i = 0; int k = 0; // for (int i = 0; i < n; i++) { while (i < matrixRow && k < matrixCol) { //看Ab[i][k]位置是否可以是主元 int maxRow = maxRow(i, k, matrixRow, addr); T[] temp = (T[]) new Number[Ab.colNum()]; //每一行无论是否有主元，都与主元最大值所在行进行交换 for (int j = 0; j < Ab.colNum(); j++) { temp[j] = Ab.get(new int[]{i, j}); Ab.set(new int[]{i, j}, Ab.get(new int[]{maxRow, j})); Ab.set(new int[]{maxRow, j}, temp[j]); } if (addr.equals(Ab.get(new int[]{i,k}),addr.zero())) { k++; }else { T u = Ab.get(new int[]{i, k}); //将主元归为一 for (int j = 0; j < Ab.colNum(); j++) { Ab.set(new int[]{i, j}, addr.div(Ab.get(new int[]{i, j}), u)); } //将主元所在行以下的行的所有非主元数据全部变成0 for (int m = i + 1; m < matrixRow; m++) { T u1 = Ab.get(new int[]{m, k}); for (int l = 0; l < Ab.colNum(); l++) { Ab.set(new int[]{m, l}, addr.sub(Ab.get(new int[]{m, l}), addr.mul(u1, Ab.get(new int[]{i, l})))); } } pivots.add(k); i++; } } } /** * 后向 */ private void backword(Addr<T> addr) { int n = pivots.size(); for (int i = n - 1; i >= 0; i--) { int k = pivots.get(i); //Ab[i][k]为主元 for (int j = i - 1; j >= 0; j--) { T u = Ab.get(new int[]{j, k}); for (int l = 0; l < Ab.colNum(); l++) { Ab.set(new int[]{j,l},addr.sub(Ab.get(new int[]{j,l}), addr.mul(u,Ab.get(new int[]{i,l})))); } } } } /** * 获取主元最大值所在的行 * @param addr * @return */ private int maxRow(int indexI,int indexJ,int n,Addr<T> addr) { T best = Ab.get(new int[]{indexI,indexJ}); int ret = indexI; for (int i = indexI + 1; i < n; i++) { if (addr.compair(Ab.get(new int[]{i,indexJ}),best) == 1) { best = Ab.get(new int[]{i,indexJ}); ret = i; } } return ret; } }
矩阵实现类
public class Matrix<T extends Number> implements Ranks<T>,Cloneable { private T[][] values; public Matrix(T[][] values) { this.values = values; } @Override public String toString() { StringBuilder builder = new StringBuilder(); for (int i = 0; i < values.length - 1; i++) { builder.append(Arrays.toString(values[i]) + ","); } builder.append(Arrays.toString(values[values.length - 1])); return "Matrix{" + "values=" + "[" + builder.toString() + "]" + '}'; } /** * 创建一个零矩阵 * @param row * @param col * @param addr * @param <T> * @return */ @SuppressWarnings("unchecked") public static <T extends Number> Ranks<T> zero(int row,int col,Addr<T> addr) { T[][] values = (T[][])new Number[row][col]; for (int i = 0; i < row; i++) { Arrays.fill(values[i],addr.zero()); } return new Matrix<>(values); } /** * 创建一个单位矩阵 * @param n * @param <T> * @return */ @SuppressWarnings("unchecked") public static <T extends Number> Ranks<T> identity(int n,Addr<T> addr) { T[][] values = (T[][])new Number[n][n]; for (int i = 0; i < n; i++) { Arrays.fill(values[i],addr.zero()); } for (int i = 0; i < n; i++) { values[i][i] = addr.one(); } return new Matrix<>(values); } @Override public int[] shape() { int[] res = new int[2]; res[0] = values.length; res[1] = values[0].length; return res; } @Override public int rowNum() { return shape()[0]; } @Override public int colNum() { return shape()[1]; } @Override public int size() { int row = rowNum(); int col = colNum(); return row * col; } @Override public T get(int pos) { if (pos >= size() || pos < 0) { throw new IllegalArgumentException("位置超出索引范围"); } int row = pos / colNum(); int col = pos % colNum(); return values[row][col]; } @Override public T get(int[] pos) { if (pos.length != 2) { throw new IllegalArgumentException("位置长度超出范围"); } if (pos[0] >= rowNum() || pos[0] < 0 || pos[1] >= colNum() || pos[1] < 0) { throw new IllegalArgumentException("位置超出索引范围"); } return values[pos[0]][pos[1]]; } @Override public void set(int[] pos, T value) { if (pos.length != 2) { throw new IllegalArgumentException("位置长度超出范围"); } if (pos[0] >= rowNum() || pos[0] < 0 || pos[1] >= colNum() || pos[1] < 0) { throw new IllegalArgumentException("位置超出索引范围"); } values[pos[0]][pos[1]] = value; } @Override public Attribute<T> rowVector(int index) { if (index >= rowNum() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } return new Vector<>(values[index]); } @Override @SuppressWarnings("unchecked") public Attribute<T> colVector(int index) { if (index >= colNum() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } T[] values = (T[])new Number[rowNum()]; for (int i = 0; i < rowNum(); i++) { values[i] = this.values[i][index]; } return new Vector<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> add(Ranks<T> another,Addr<T> addr) { if (!Arrays.equals(this.shape(),another.shape())) { throw new IllegalArgumentException("加法错误，两个矩阵的形状必须相同"); } T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.add(this.get(new int[]{i,j}),another.get(new int[]{i,j})); } } this.values = values; return this; } @Override @SuppressWarnings("unchecked") public Ranks<T> sub(Ranks<T> another, Addr<T> addr) { if (!Arrays.equals(this.shape(),another.shape())) { throw new IllegalArgumentException("减法错误，两个矩阵的形状必须相同"); } T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.sub(this.get(new int[]{i,j}),another.get(new int[]{i,j})); } } this.values = values; return this; } @Override @SuppressWarnings("unchecked") public Ranks<T> mul(T k, Addr<T> addr) { T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.mul(k,this.get(new int[]{i,j})); } } this.values = values; return this; } @Override public Ranks<Double> div(double k, Addr<T> addr) { Double[][] values = new Double[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.div(this.get(new int[]{i,j}),k); } } return new Matrix<>(values); } @Override public Ranks<T> pos(Addr<T> addr,T one) { return mul(one,addr); } @Override public Ranks<T> neg(Addr<T> addr,T negOne) { return mul(negOne,addr); } @Override @SuppressWarnings("unchecked") public Attribute<T> dot(Attribute<T> another, Addr<T> addr) { if (this.colNum() != another.len()) { throw new IllegalArgumentException("矩阵列数需要与向量的长度相等"); } T[] values = (T[])new Number[this.rowNum()]; for (int i = 0; i < this.rowNum(); i++) { values[i] = this.rowVector(i).dot(another,addr); } return new Vector<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> dot(Ranks<T> another, Addr<T> addr) { if (this.colNum() != another.rowNum()) { throw new IllegalArgumentException("当前矩阵列数必须等于点乘矩阵行数"); } T[][] values = (T[][]) new Number[this.rowNum()][another.colNum()]; for (int i = 0; i < this.rowNum(); i++) { for (int j = 0; j < another.colNum(); j++) { values[i][j] = this.rowVector(i).dot(another.colVector(j),addr); } } return new Matrix<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> trans() { T[][] values = (T[][])new Number[this.colNum()][this.rowNum()]; for (int i = 0; i < this.colNum(); i++) { for (int j = 0; j < this.rowNum(); j++) { values[i][j] = this.values[j][i]; } } return new Matrix<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> inv(Addr<T> addr) { if (this.rowNum() != this.colNum()) { return null; } int n = this.rowNum(); FuncGroup<T> ls = new LinearSystem<>(this,Matrix.identity(n,addr)); if (!ls.gaussJordanElimination(addr)) { return null; } T[][] values = (T[][])new Number[this.rowNum()][((LinearSystem)ls).getAb().colNum() - this.colNum()]; for (int i = 0; i < this.rowNum(); i++) { for (int j = 0; j < ((LinearSystem)ls).getAb().colNum() - this.colNum(); j++) { values[i][j] = (T)((LinearSystem)ls).getAb().get(new int[]{i,((LinearSystem)ls).getAb().colNum() - this.colNum() + j}); } } return new Matrix<>(values); } @Override @SuppressWarnings("unchecked") public Matrix clone() throws CloneNotSupportedException { T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = this.values[i][j]; } } return new Matrix(values); } @SuppressWarnings("unchecked") public static void main(String[] args) throws CloneNotSupportedException { Ranks<Integer> matrix = new Matrix<>(new Integer[][]{{1,2,3},{3,4,5}}); System.out.println(matrix); System.out.println(Arrays.toString(matrix.shape())); System.out.println(matrix.rowNum()); System.out.println(matrix.colNum()); System.out.println(matrix.size()); System.out.println(matrix.get(4)); System.out.println(matrix.get(new int[]{0,0})); System.out.println(matrix.rowVector(1)); System.out.println(matrix.colVector(2)); Ranks<Integer> matrix1 = new Matrix<>(new Integer[][]{{7,8,9},{10,20,30}}); Addr<Integer> addr = new Addr<Integer>() { @Override public Integer zero() { return 0; } @Override public Integer one() { return 1; } @Override public Integer add(Integer lhs, Integer rhs) { return lhs + rhs; } @Override public Integer sub(Integer lhs, Integer rhs) { return lhs - rhs; } @Override public Integer mul(Integer k, Integer hs) { return k * hs; } @Override public Integer square(Integer hs) { return hs * hs; } @Override public double prescription(Integer hs) { return Math.sqrt(hs); } @Override public double div(Integer hs, double nhs) { return hs / nhs; } @Override public Integer div(Integer hs, Integer nhs) { return hs / nhs; } @Override public double acos(double hs) { return Math.acos(hs); } @Override public double toDegree(double hs) { return Math.toDegrees(hs); } @Override public int compair(Integer lhs, Integer rhs) { if (lhs - rhs > 0) { return 1; }else if (lhs - rhs < 0) { return -1; }else { return 0; } } @Override public boolean equals(Integer lhs, Integer rhs) { return lhs == rhs; } }; System.out.println(matrix.trans()); Ranks<Integer> copy = ((Matrix)matrix).clone(); Attribute<Integer> vector = new Vector<>(new Integer[]{9,33,97}); Ranks<Integer> matrix2 = new Matrix<>(new Integer[][]{{7,9},{11,33},{42,97}}); //[(1,2,3),(3,4,5)]*(23,65,79) System.out.println(matrix.dot(vector,addr)); // [(1,2,3),(3,4,5)]*[(7,9),(11,33),(42,97)] System.out.println(matrix.dot(matrix2,addr)); Ranks<Integer> copy1 = ((Matrix)matrix).clone(); System.out.println(matrix.add(matrix1,addr)); System.out.println(matrix1.sub(copy,addr)); System.out.println(copy.div(2,addr)); System.out.println(copy.pos(addr,1)); System.out.println(copy.neg(addr,-1)); System.out.println(copy1.mul(3,addr)); System.out.println(Matrix.zero(2,3,addr)); Ranks<Integer> I = Matrix.identity(4, addr); System.out.println(I); Matrix<Integer> idenTest = new Matrix<>(new Integer[][]{{1,2,3,4},{5,6,7,8},{9,10,11,12},{13,14,15,16}}); System.out.println(idenTest.dot(I,addr)); System.out.println(I.dot(idenTest,addr)); Addr<Double> addr1 = new Addr<Double>() { @Override public Double zero() { return 0.0; } @Override public Double one() { return 1.0; } @Override public Double add(Double lhs, Double rhs) { return lhs + rhs; } @Override public Double sub(Double lhs, Double rhs) { return lhs - rhs; } @Override public Double mul(Double k, Double hs) { return k * hs; } @Override public Double square(Double hs) { return hs * hs; } @Override public double prescription(Double hs) { return Math.sqrt(hs); } @Override public double div(Double hs, double nhs) { return hs / nhs; } @Override public Double div(Double hs, Double nhs) { return hs / nhs; } @Override public double acos(double hs) { return Math.acos(hs); } @Override public double toDegree(double hs) { return Math.toDegrees(hs); } @Override public int compair(Double lhs, Double rhs) { if (lhs - rhs > 0) { return 1; }else if (lhs - rhs < 0) { return -1; }else { return 0; } } @Override public boolean equals(Double lhs, Double rhs) { double dis = 1e-8; if (Math.abs(lhs - rhs) < dis) { return true; } return false; } }; Ranks<Double> A = new Matrix<>(new Double[][]{{1.0,2.0},{3.0,4.0}}); System.out.println(A.inv(addr1)); } }
运行结果
Matrix{values=[[1, 2, 3],[3, 4, 5]]}
[2, 3]
2
3
6
4
1
Vector{values=[3, 4, 5]}
Vector{values=[3, 5]}
Matrix{values=[[1, 3],[2, 4],[3, 5]]}
Vector{values=[366, 644]}
Matrix{values=[[155, 366],[275, 644]]}
Matrix{values=[[8, 10, 12],[13, 24, 35]]}
Matrix{values=[[6, 6, 6],[7, 16, 25]]}
Matrix{values=[[0.5, 1.0, 1.5],[1.5, 2.0, 2.5]]}
Matrix{values=[[1, 2, 3],[3, 4, 5]]}
Matrix{values=[[-1, -2, -3],[-3, -4, -5]]}
Matrix{values=[[3, 6, 9],[9, 12, 15]]}
Matrix{values=[[0, 0, 0],[0, 0, 0]]}
Matrix{values=[[1, 0, 0, 0],[0, 1, 0, 0],[0, 0, 1, 0],[0, 0, 0, 1]]}
Matrix{values=[[1, 2, 3, 4],[5, 6, 7, 8],[9, 10, 11, 12],[13, 14, 15, 16]]}
Matrix{values=[[1, 2, 3, 4],[5, 6, 7, 8],[9, 10, 11, 12],[13, 14, 15, 16]]}
Matrix{values=[[-1.9999999999999996, 0.9999999999999998],[1.4999999999999998, -0.49999999999999994]]}
Python代码
from .Vector import Vector from .Matrix import Matrix from ._global import is_zero class LinearSystem: def __init__(self, A, b): assert A.row_num() == len(b), "矩阵的行数必须等于向量的长度" self._matrix_row = A.row_num() self._matrix_col = A.col_num() if isinstance(b, Vector): self.Ab = [Vector(A.row_vector(i).underlying_list() + [b[i]]) for i in range(self._matrix_row)] if isinstance(b, Matrix): self.Ab = [Vector(A.row_vector(i).underlying_list() + b.row_vector(i).underlying_list()) for i in range(self._matrix_row)] self.pivots = [] def _max_row(self, index_i, index_j, n): # 主元最大值所在的行 best, ret = self.Ab[index_i][index_j], index_i for i in range(index_i + 1, n): if self.Ab[i][index_j] > best: best, ret = self.Ab[i][index_j], i return ret def _forward(self): # 前向 i, k = 0, 0 while i < self._matrix_row and k < self._matrix_col: # 看Ab[i][k]位置是否可以为主元 max_row = self._max_row(i, k, self._matrix_row) self.Ab[i], self.Ab[max_row] = self.Ab[max_row], self.Ab[i] if is_zero(self.Ab[i][k]): k += 1 else: # 将主元归为一 self.Ab[i] = self.Ab[i] / self.Ab[i][k] for j in range(i + 1, self._matrix_row): self.Ab[j] = self.Ab[j] - self.Ab[j][k] * self.Ab[i] self.pivots.append(k) i += 1 def _backward(self): # 后向 n = len(self.pivots) for i in range(n - 1, -1, -1): k = self.pivots[i] # Ab[i][k]为主元 for j in range(i - 1, -1, -1): self.Ab[j] = self.Ab[j] - self.Ab[j][k] * self.Ab[i] def gauss_jordan_elimination(self): # 高斯-约旦消元法，如果有解，返回True;如果没有解，返回False self._forward() self._backward() for i in range(len(self.pivots), self._matrix_row): if not is_zero(self.Ab[i][-1]): return False return True def fancy_print(self): # 打印增广矩阵 for i in range(self._matrix_row): print(" ".join(str(self.Ab[i][j]) for j in range(self._matrix_col)), end=" ") print("|", self.Ab[i][-1])
from .Vector import Vector class Matrix: def __init__(self, list2d): self._values = [row[:] for row in list2d] @classmethod def zero(cls, r, c): # 返回一个r行c列的零矩阵 return cls([[0] * c for _ in range(r)]) @classmethod def identity(cls, n): # 返回一个n行n列的单位矩阵 m = [[0] * n for _ in range(n)] for i in range(n): m[i][i] = 1 return cls(m) def trans(self): # 返回矩阵的转置矩阵 return Matrix([[e for e in self.col_vector(i)] for i in range(self.col_num())]) def inv(self): # 返回矩阵的逆 from .LinearSystem import LinearSystem if self.row_num() != self.col_num(): return None n = self.row_num() ls = LinearSystem(self, Matrix.identity(n)) if not ls.gauss_jordan_elimination(): return None invA = [[row[i] for i in range(n, 2 * n)] for row in ls.Ab] return Matrix(invA) def __add__(self, another): # 返回两个矩阵的加法结果 assert self.shape() == another.shape(), \ "加法错误，两个矩阵的形状必须相同" return Matrix([[a + b for a, b in zip(self.row_vector(i), another.row_vector(i))] for i in range(self.row_num())]) def __sub__(self, another): # 返回两个矩阵的减法结果 assert self.shape() == another.shape(), \ "减法错误，两个矩阵的形状必须相同" return Matrix([[a - b for a, b in zip(self.row_vector(i), another.row_vector(i))] for i in range(self.row_num())]) def dot(self, another): if isinstance(another, Vector): # 矩阵和向量的乘法 assert self.col_num() == len(another), \ "矩阵列数需要与向量的长度相等" return Vector([self.row_vector(i).dot(another) for i in range(self.row_num())]) if isinstance(another, Matrix): # 矩阵和矩阵的乘法 assert self.col_num() == another.row_num(), \ "当前矩阵列数必须等于点乘矩阵行数" return Matrix([[self.row_vector(i).dot(another.col_vector(j)) for j in range(another.col_num())] for i in range(self.row_num())]) def __mul__(self, k): # 返回矩阵的数量乘结果: self * k return Matrix([[e * k for e in self.row_vector(i)] for i in range(self.row_num())]) def __rmul__(self, k): # 返回矩阵的数量乘结果: k * self return self * k def __truediv__(self, k): # 返回数量除法的结果矩阵: self / k return (1 / k) * self def __pos__(self): # 返回矩阵取正的结果 return 1 * self def __neg__(self): # 返回矩阵取负的结果 return -1 * self def row_vector(self, index): # 返回矩阵的第index个行向量 return Vector(self._values[index]) def col_vector(self, index): # 返回矩阵的第index个列向量 return Vector([row[index] for row in self._values]) def __getitem__(self, pos): # 返回矩阵pos位置的元素 r, c = pos return self._values[r][c] def size(self): # 返回矩阵中元素的个数 r, c = self.shape() return r * c def row_num(self): # 返回矩阵的行数 return self.shape()[0] __len__ = row_num def col_num(self): # 返回矩阵的列数 return self.shape()[1] def shape(self): # 返回矩阵的形状：(行数，列数) return len(self._values), len(self._values[0]) def __repr__(self): return "Matrix({})".format(self._values) __str__ = __repr__
from playLA.Matrix import Matrix from playLA.Vector import Vector if __name__ == "__main__": matrix = Matrix([[1, 2, 3], [3, 4, 5]]) print(matrix) print("matrix.shape = {}".format(matrix.shape())) print("matrix.size = {}".format(matrix.size())) print("len(matrix) = {}".format(len(matrix))) print("matrix[0][0] = {}".format(matrix[0, 0])) print(matrix.row_vector(1)) print(matrix.col_vector(0)) matrix2 = Matrix([[7, 8, 9], [10, 20, 30]]) print("add: {}".format(matrix + matrix2)) print("subtract: {}".format(matrix2 - matrix)) print("scalar-mul: {}".format(matrix * 3)) print("scalar-mul: {}".format(3 * matrix)) print("zero_2_3: {}".format(Matrix.zero(2, 3))) vec = Vector([23, 65, 79]) print("matrix.dot(vec) = {}".format(matrix.dot(vec))) matrix3 = Matrix([[7, 9], [11, 33], [42, 97]]) print("matrix.dot(matrix3) = {}".format(matrix.dot(matrix3))) print("matrix.trans = {}".format(matrix.trans())) I = Matrix.identity(4) print(I) idenTest = Matrix([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]]) print(idenTest.dot(I)) print(I.dot(idenTest)) A = Matrix([[1, 2], [3, 4]]) print(A.inv())
运行结果
Matrix([[1, 2, 3], [3, 4, 5]])
matrix.shape = (2, 3)
matrix.size = 6
len(matrix) = 2
matrix[0][0] = 1
(3, 4, 5)
(1, 3)
add: Matrix([[8, 10, 12], [13, 24, 35]])
subtract: Matrix([[6, 6, 6], [7, 16, 25]])
scalar-mul: Matrix([[3, 6, 9], [9, 12, 15]])
scalar-mul: Matrix([[3, 6, 9], [9, 12, 15]])
zero_2_3: Matrix([[0, 0, 0], [0, 0, 0]])
matrix.dot(vec) = (390, 724)
matrix.dot(matrix3) = Matrix([[155, 366], [275, 644]])
matrix.trans = Matrix([[1, 3], [2, 4], [3, 5]])
Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])
Matrix([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]])
Matrix([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]])
Matrix([[-1.9999999999999996, 0.9999999999999998], [1.4999999999999998, -0.4999999999999999]])
初等矩阵
在我们使用高斯-约旦消元法中，我们会矩阵的某一行乘以一个常数，矩阵的一行加(减)零一行，交换矩阵的两行
这些都是在矩阵的内部，对内部元素进行操作。但矩阵本身表示变换，对于这些操作，同样可以使用矩阵来表示
我们知道单位矩阵乘以一个矩阵等于该矩阵本身，其意义就是单位矩阵作为一个变换矩阵，没有对该矩阵产生变换而已。
如果想让矩阵的某一行乘以一个常数k，我们可以对单位矩阵做出如下的改变(比如说让矩阵的第一行乘以k)
如果想让矩阵的一行加(减)另一行，我们可以对单位矩阵做出如下改变(比如说让矩阵的第三行加上第一行)
再比如说让矩阵的第一行减去第三行，我们可以对单位矩阵做出如下改变
当然我们也可以加(减)另一行的若干倍,我们可以对单位矩阵做出如下改变(比如说让矩阵的第二行减去p倍的第三行)
现在我们想交换矩阵的两行，我们可以对单位矩阵做出如下改变(比如说交换矩阵的第二行和第三行)
以上的这些变换过程，是矩阵的基本操作，称为矩阵的初等变换，而这些单位矩阵变换后的矩阵，称为初等矩阵
则初等矩阵就是形如下面这样的矩阵
而初等矩阵就是对单位矩阵进行一次初等变换得到的结果矩阵，通常记作E
我们在高斯-约旦消元法把矩阵化为行最简形式的过程就是对增广矩阵A进行了一系列的初等变换，就相当于是寻找一系列的初等矩阵E，使得：
其中rref(A)表示A的行最简形式。
Addr接口增加一个取负的方法neg()
public interface Addr<T extends Number> { T zero(); T one(); T add(T lhs,T rhs); T sub(T lhs,T rhs); T mul(T k,T hs); T square(T hs); double prescription(T hs); double div(T hs,double nhs); T div(T hs,T nhs); double acos(double hs); double toDegree(double hs); int compair(T lhs,T rhs); boolean equals(T lhs,T rhs); T neg(T hs); }
我们将之前的线性系统类改成一个抽象类
public abstract class LinearSystem<T extends Number> implements FuncGroup<T> { //系数矩阵的行数 protected int matrixRow; //系数矩阵的列数 protected int matrixCol; //增广矩阵 @Getter protected Ranks<T> Ab; //主元列 @Getter protected List<Integer> pivots = new ArrayList<>(); @SuppressWarnings("unchecked") public LinearSystem(Ranks<T> A, Attribute<T> b) { if (A.rowNum() != b.len()) { throw new IllegalArgumentException("系数矩阵的行数必须等于常数向量的长度"); } matrixRow = A.rowNum(); matrixCol = A.colNum(); T[][] values = (T[][])new Number[A.rowNum()][A.colNum() + 1]; for (int i = 0; i < A.rowNum(); i++) { for (int j = 0; j <= A.colNum(); j++) { if (j < A.colNum()) { values[i][j] = A.get(new int[]{i, j}); }else { values[i][j] = b.get(i); } } } Ab = new Matrix<>(values); } @SuppressWarnings("unchecked") public LinearSystem(Ranks<T> A, Ranks<T> b) { if (A.rowNum() != b.colNum()) { throw new IllegalArgumentException("系数矩阵的行数必须等于单位矩阵的列数"); } matrixRow = A.rowNum(); matrixCol = A.colNum(); T[][] values = (T[][])new Number[A.rowNum()][A.colNum() + b.colNum()]; for (int i = 0; i < A.rowNum(); i++) { for (int j = 0; j < A.colNum() + b.colNum(); j++) { if (j < A.colNum()) { values[i][j] = A.get(new int[]{i, j}); }else { values[i][j] = b.get(new int[]{i, j - A.colNum()}); } } } Ab = new Matrix<>(values); } @Override public boolean gaussJordanElimination(Addr<T> addr) { forward(addr); backword(addr); //无主元的行的非系数元素如果为0则无解 for (int i = pivots.size(); i < matrixRow; i++) { if (!addr.equals(Ab.get(new int[]{i,matrixCol}),addr.zero())) { return false; } } return true; } @Override public void fancyPrint() { StringBuilder builder = new StringBuilder(); if (Ab.colNum() - matrixCol > 1) { for (int i = 0; i < matrixRow; i++) { for (int j = 0; j < Ab.colNum(); j++) { if (j < matrixCol) { builder.append(" "); builder.append(Ab.get(new int[]{i, j})); } else if (j == matrixCol) { builder.append(" | "); builder.append(Ab.get(new int[]{i, j})); } else { builder.append(" "); builder.append(Ab.get(new int[]{i, j})); builder.append("\n"); } } } }else if (Ab.colNum() - matrixCol == 1) { for (int i = 0; i < matrixRow; i++) { for (int j = 0; j < Ab.colNum(); j++) { if (j < matrixCol) { builder.append(" "); builder.append(Ab.get(new int[]{i, j})); } else if (j == matrixCol) { builder.append(" | "); builder.append(Ab.get(new int[]{i, j})); builder.append("\n"); } } } } System.out.println(builder.toString()); } protected abstract void forward(Addr<T> addr); protected abstract void backword(Addr<T> addr); /** * 获取主元最大值所在的行 * @param addr * @return */ protected int maxRow(int indexI,int indexJ,int n,Addr<T> addr) { T best = Ab.get(new int[]{indexI,indexJ}); int ret = indexI; for (int i = indexI + 1; i < n; i++) { if (addr.compair(Ab.get(new int[]{i,indexJ}),best) == 1) { best = Ab.get(new int[]{i,indexJ}); ret = i; } } return ret; } }
之前的增广矩阵内部元素的处理类改成LinearSystemFunc类
/** * 线性系统 * @param <T> */ public final class LinearSystemFunc<T extends Number> extends LinearSystem<T>{ public LinearSystemFunc(Ranks<T> A, Attribute<T> b) { super(A,b); } public LinearSystemFunc(Ranks<T> A, Ranks<T> b) { super(A,b); } /** * 前向 * @param addr */ @Override @SuppressWarnings("unchecked") protected void forward(Addr<T> addr) { // int n = matrixRow; int i = 0; int k = 0; // for (int i = 0; i < n; i++) { while (i < matrixRow && k < matrixCol) { //看Ab[i][k]位置是否可以是主元 int maxRow = maxRow(i, k, matrixRow, addr); T[] temp = (T[]) new Number[Ab.colNum()]; //每一行无论是否有主元，都与主元最大值所在行进行交换 for (int j = 0; j < Ab.colNum(); j++) { temp[j] = Ab.get(new int[]{i, j}); Ab.set(new int[]{i, j}, Ab.get(new int[]{maxRow, j})); Ab.set(new int[]{maxRow, j}, temp[j]); } if (addr.equals(Ab.get(new int[]{i,k}),addr.zero())) { k++; }else { T u = Ab.get(new int[]{i, k}); //将主元归为一 for (int j = 0; j < Ab.colNum(); j++) { Ab.set(new int[]{i, j}, addr.div(Ab.get(new int[]{i, j}), u)); } //将主元所在行以下的行的所有非主元数据全部变成0 for (int m = i + 1; m < matrixRow; m++) { T u1 = Ab.get(new int[]{m, k}); for (int l = 0; l < Ab.colNum(); l++) { Ab.set(new int[]{m, l}, addr.sub(Ab.get(new int[]{m, l}), addr.mul(u1, Ab.get(new int[]{i, l})))); } } pivots.add(k); i++; } } } /** * 后向 */ @Override protected void backword(Addr<T> addr) { int n = pivots.size(); for (int i = n - 1; i >= 0; i--) { int k = pivots.get(i); //Ab[i][k]为主元 for (int j = i - 1; j >= 0; j--) { T u = Ab.get(new int[]{j, k}); for (int l = 0; l < Ab.colNum(); l++) { Ab.set(new int[]{j,l},addr.sub(Ab.get(new int[]{j,l}), addr.mul(u,Ab.get(new int[]{i,l})))); } } } } }
新增加我们的初等矩阵变换的线形系统类LinearSystemMatrix
/** * 初等矩阵变换线性系统 * @param <T> */ public final class LinearSystemMatrix<T extends Number> extends LinearSystem<T> { public LinearSystemMatrix(Ranks<T> A,Attribute<T> b) { super(A,b); } public LinearSystemMatrix(Ranks<T> A,Ranks<T> b) { super(A,b); } @Override protected void forward(Addr<T> addr) { int i = 0; int k = 0; while (i < matrixRow && k < matrixCol) { int maxRow = maxRow(i, k, matrixRow, addr); Ab = changeMatrix(i,maxRow,addr).dot(Ab,addr); if (addr.equals(Ab.get(new int[]{i,k}),addr.zero())) { k++; }else { T u = Ab.get(new int[]{i, k}); Ab = mulMatrix(i,addr.div(addr.one(),u),addr).dot(Ab,addr); for (int j = i + 1; j < matrixRow; j++) { T u1 = Ab.get(new int[]{j, k}); Ab = addPMatrix(j,i,addr.neg(u1),addr).dot(Ab,addr); } pivots.add(k); i++; } } } @Override protected void backword(Addr<T> addr) { int n = pivots.size(); for (int i = n - 1; i >= 0; i--) { int k = pivots.get(i); //Ab[i][k]为主元 for (int j = i - 1; j >= 0; j--) { T u = Ab.get(new int[]{j, k}); Ab = addPMatrix(j,i,addr.neg(u),addr).dot(Ab,addr); } } } /** * 交换初等矩阵 * @param i 原行 * @param j 交换行 * @param addr * @return */ private Ranks<T> changeMatrix(int i,int j,Addr<T> addr) { Ranks<T> identity = Matrix.identity(matrixRow,addr); identity.set(new int[]{i,i},addr.zero()); identity.set(new int[]{j,i},addr.one()); identity.set(new int[]{j,j},addr.zero()); identity.set(new int[]{i,j},addr.one()); return identity; } /** * 乘除初等矩阵 * @param i 需要处理的矩阵的行 * @param t 乘的系数 * @param addr * @return */ private Ranks<T> mulMatrix(int i,T t,Addr<T> addr) { Ranks<T> identity = Matrix.identity(matrixRow,addr); identity.set(new int[]{i,i},t); return identity; } /** * 加减初等矩阵 * @param i 需要加减的行 * @param j 被加减的行 * @param t 加减的倍数 * @param addr * @return */ private Ranks<T> addPMatrix(int i,int j,T t,Addr<T> addr) { Ranks<T> identity = Matrix.identity(matrixRow,addr); identity.set(new int[]{i,j},t); return identity; } }
将Ranks接口的求矩阵的逆改成两个方法
public interface Ranks<T extends Number> { /** * 矩阵的行列数 * @return */ int[] shape(); /** * 矩阵的行数 * @return */ int rowNum(); /** * 矩阵的列数 * @return */ int colNum(); /** * 矩阵的元素个数 * @return */ int size(); /** * 获取矩阵中的一个元素 * @param pos * @return */ T get(int pos); /** * 获取矩阵中的一个元素 * @param pos * @return */ T get(int[] pos); /** * 设置矩阵中某个位置的元素 * @param pos * @param value */ void set(int[] pos,T value); /** * 获取一个行向量 * @param index * @return */ Attribute<T> rowVector(int index); /** * 获取一个列向量 * @param index * @return */ Attribute<T> colVector(int index); /** * 矩阵相加 * @param another * @return */ Ranks<T> add(Ranks<T> another,Addr<T> addr); /** * 矩阵相减 * @param another * @param addr * @return */ Ranks<T> sub(Ranks<T> another,Addr<T> addr); /** * 矩阵乘以一个数 * @param k * @param addr * @return */ Ranks<T> mul(T k,Addr<T> addr); /** * 矩阵除以一个数 * @param k * @param addr * @return */ Ranks<Double> div(double k,Addr<T> addr); /** * 矩阵取正的结果 * @return */ Ranks<T> pos(Addr<T> addr,T one); /** * 矩阵取负的结果 * @return */ Ranks<T> neg(Addr<T> addr,T negOne); /** * 矩阵与向量的点乘 * @param another * @param addr * @return */ Attribute<T> dot(Attribute<T> another,Addr<T> addr); /** * 矩阵与矩阵的点乘 * @param another * @param addr * @return */ Ranks<T> dot(Ranks<T> another,Addr<T> addr); /** * 转置矩阵 * @return */ Ranks<T> trans(); /** * 矩阵的逆 * @param addr * @return */ Ranks<T> invFunc(Addr<T> addr); /** * 矩阵的逆 * @param addr * @return */ Ranks<T> invMatrix(Addr<T> addr); }
矩阵实现类
public class Matrix<T extends Number> implements Ranks<T>,Cloneable { private T[][] values; public Matrix(T[][] values) { this.values = values; } @Override public String toString() { StringBuilder builder = new StringBuilder(); for (int i = 0; i < values.length - 1; i++) { builder.append(Arrays.toString(values[i]) + ","); } builder.append(Arrays.toString(values[values.length - 1])); return "Matrix{" + "values=" + "[" + builder.toString() + "]" + '}'; } /** * 创建一个零矩阵 * @param row * @param col * @param addr * @param <T> * @return */ @SuppressWarnings("unchecked") public static <T extends Number> Ranks<T> zero(int row,int col,Addr<T> addr) { T[][] values = (T[][])new Number[row][col]; for (int i = 0; i < row; i++) { Arrays.fill(values[i],addr.zero()); } return new Matrix<>(values); } /** * 创建一个单位矩阵 * @param n * @param <T> * @return */ @SuppressWarnings("unchecked") public static <T extends Number> Ranks<T> identity(int n,Addr<T> addr) { T[][] values = (T[][])new Number[n][n]; for (int i = 0; i < n; i++) { Arrays.fill(values[i],addr.zero()); } for (int i = 0; i < n; i++) { values[i][i] = addr.one(); } return new Matrix<>(values); } @Override public int[] shape() { int[] res = new int[2]; res[0] = values.length; res[1] = values[0].length; return res; } @Override public int rowNum() { return shape()[0]; } @Override public int colNum() { return shape()[1]; } @Override public int size() { int row = rowNum(); int col = colNum(); return row * col; } @Override public T get(int pos) { if (pos >= size() || pos < 0) { throw new IllegalArgumentException("位置超出索引范围"); } int row = pos / colNum(); int col = pos % colNum(); return values[row][col]; } @Override public T get(int[] pos) { if (pos.length != 2) { throw new IllegalArgumentException("位置长度超出范围"); } if (pos[0] >= rowNum() || pos[0] < 0 || pos[1] >= colNum() || pos[1] < 0) { throw new IllegalArgumentException("位置超出索引范围"); } return values[pos[0]][pos[1]]; } @Override public void set(int[] pos, T value) { if (pos.length != 2) { throw new IllegalArgumentException("位置长度超出范围"); } if (pos[0] >= rowNum() || pos[0] < 0 || pos[1] >= colNum() || pos[1] < 0) { throw new IllegalArgumentException("位置超出索引范围"); } values[pos[0]][pos[1]] = value; } @Override public Attribute<T> rowVector(int index) { if (index >= rowNum() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } return new Vector<>(values[index]); } @Override @SuppressWarnings("unchecked") public Attribute<T> colVector(int index) { if (index >= colNum() || index < 0) { throw new IllegalArgumentException("索引超出范围"); } T[] values = (T[])new Number[rowNum()]; for (int i = 0; i < rowNum(); i++) { values[i] = this.values[i][index]; } return new Vector<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> add(Ranks<T> another,Addr<T> addr) { if (!Arrays.equals(this.shape(),another.shape())) { throw new IllegalArgumentException("加法错误，两个矩阵的形状必须相同"); } T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.add(this.get(new int[]{i,j}),another.get(new int[]{i,j})); } } this.values = values; return this; } @Override @SuppressWarnings("unchecked") public Ranks<T> sub(Ranks<T> another, Addr<T> addr) { if (!Arrays.equals(this.shape(),another.shape())) { throw new IllegalArgumentException("减法错误，两个矩阵的形状必须相同"); } T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.sub(this.get(new int[]{i,j}),another.get(new int[]{i,j})); } } this.values = values; return this; } @Override @SuppressWarnings("unchecked") public Ranks<T> mul(T k, Addr<T> addr) { T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.mul(k,this.get(new int[]{i,j})); } } this.values = values; return this; } @Override public Ranks<Double> div(double k, Addr<T> addr) { Double[][] values = new Double[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = addr.div(this.get(new int[]{i,j}),k); } } return new Matrix<>(values); } @Override public Ranks<T> pos(Addr<T> addr,T one) { return mul(one,addr); } @Override public Ranks<T> neg(Addr<T> addr,T negOne) { return mul(negOne,addr); } @Override @SuppressWarnings("unchecked") public Attribute<T> dot(Attribute<T> another, Addr<T> addr) { if (this.colNum() != another.len()) { throw new IllegalArgumentException("矩阵列数需要与向量的长度相等"); } T[] values = (T[])new Number[this.rowNum()]; for (int i = 0; i < this.rowNum(); i++) { values[i] = this.rowVector(i).dot(another,addr); } return new Vector<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> dot(Ranks<T> another, Addr<T> addr) { if (this.colNum() != another.rowNum()) { throw new IllegalArgumentException("当前矩阵列数必须等于点乘矩阵行数"); } T[][] values = (T[][]) new Number[this.rowNum()][another.colNum()]; for (int i = 0; i < this.rowNum(); i++) { for (int j = 0; j < another.colNum(); j++) { values[i][j] = this.rowVector(i).dot(another.colVector(j),addr); } } return new Matrix<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> trans() { T[][] values = (T[][])new Number[this.colNum()][this.rowNum()]; for (int i = 0; i < this.colNum(); i++) { for (int j = 0; j < this.rowNum(); j++) { values[i][j] = this.values[j][i]; } } return new Matrix<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> invFunc(Addr<T> addr) { if (this.rowNum() != this.colNum()) { return null; } int n = this.rowNum(); FuncGroup<T> ls = new LinearSystemFunc<>(this,Matrix.identity(n,addr)); if (!ls.gaussJordanElimination(addr)) { return null; } T[][] values = (T[][])new Number[this.rowNum()][((LinearSystemFunc)ls).getAb().colNum() - this.colNum()]; for (int i = 0; i < this.rowNum(); i++) { for (int j = 0; j < ((LinearSystemFunc)ls).getAb().colNum() - this.colNum(); j++) { values[i][j] = (T)((LinearSystemFunc)ls).getAb().get(new int[]{i,((LinearSystemFunc)ls).getAb().colNum() - this.colNum() + j}); } } return new Matrix<>(values); } @Override @SuppressWarnings("unchecked") public Ranks<T> invMatrix(Addr<T> addr) { if (this.rowNum() != this.colNum()) { return null; } int n = this.rowNum(); FuncGroup<T> ls = new LinearSystemMatrix<>(this,Matrix.identity(n,addr)); if (!ls.gaussJordanElimination(addr)) { return null; } T[][] values = (T[][])new Number[this.rowNum()][((LinearSystemMatrix)ls).getAb().colNum() - this.colNum()]; for (int i = 0; i < this.rowNum(); i++) { for (int j = 0; j < ((LinearSystemMatrix)ls).getAb().colNum() - this.colNum(); j++) { values[i][j] = (T)((LinearSystemMatrix)ls).getAb().get(new int[]{i,((LinearSystemMatrix)ls).getAb().colNum() - this.colNum() + j}); } } return new Matrix<>(values); } @Override @SuppressWarnings("unchecked") public Matrix clone() throws CloneNotSupportedException { T[][] values = (T[][])new Number[this.rowNum()][this.colNum()]; for (int i = 0; i < rowNum(); i++) { for (int j = 0; j < colNum(); j++) { values[i][j] = this.values[i][j]; } } return new Matrix(values); } @SuppressWarnings("unchecked") public static void main(String[] args) throws CloneNotSupportedException { Ranks<Integer> matrix = new Matrix<>(new Integer[][]{{1,2,3},{3,4,5}}); System.out.println(matrix); System.out.println(Arrays.toString(matrix.shape())); System.out.println(matrix.rowNum()); System.out.println(matrix.colNum()); System.out.println(matrix.size()); System.out.println(matrix.get(4)); System.out.println(matrix.get(new int[]{0,0})); System.out.println(matrix.rowVector(1)); System.out.println(matrix.colVector(2)); Ranks<Integer> matrix1 = new Matrix<>(new Integer[][]{{7,8,9},{10,20,30}}); Addr<Integer> addr = new Addr<Integer>() { @Override public Integer zero() { return 0; } @Override public Integer one() { return 1; } @Override public Integer add(Integer lhs, Integer rhs) { return lhs + rhs; } @Override public Integer sub(Integer lhs, Integer rhs) { return lhs - rhs; } @Override public Integer mul(Integer k, Integer hs) { return k * hs; } @Override public Integer square(Integer hs) { return hs * hs; } @Override public double prescription(Integer hs) { return Math.sqrt(hs); } @Override public double div(Integer hs, double nhs) { return hs / nhs; } @Override public Integer div(Integer hs, Integer nhs) { return hs / nhs; } @Override public double acos(double hs) { return Math.acos(hs); } @Override public double toDegree(double hs) { return Math.toDegrees(hs); } @Override public int compair(Integer lhs, Integer rhs) { if (lhs - rhs > 0) { return 1; }else if (lhs - rhs < 0) { return -1; }else { return 0; } } @Override public boolean equals(Integer lhs, Integer rhs) { return lhs == rhs; } @Override public Integer neg(Integer hs) { return -hs; } }; System.out.println(matrix.trans()); Ranks<Integer> copy = ((Matrix)matrix).clone(); Attribute<Integer> vector = new Vector<>(new Integer[]{9,33,97}); Ranks<Integer> matrix2 = new Matrix<>(new Integer[][]{{7,9},{11,33},{42,97}}); //[(1,2,3),(3,4,5)]*(23,65,79) System.out.println(matrix.dot(vector,addr)); // [(1,2,3),(3,4,5)]*[(7,9),(11,33),(42,97)] System.out.println(matrix.dot(matrix2,addr)); Ranks<Integer> copy1 = ((Matrix)matrix).clone(); System.out.println(matrix.add(matrix1,addr)); System.out.println(matrix1.sub(copy,addr)); System.out.println(copy.div(2,addr)); System.out.println(copy.pos(addr,1)); System.out.println(copy.neg(addr,-1)); System.out.println(copy1.mul(3,addr)); System.out.println(Matrix.zero(2,3,addr)); Ranks<Integer> I = Matrix.identity(4, addr); System.out.println(I); Matrix<Integer> idenTest = new Matrix<>(new Integer[][]{{1,2,3,4},{5,6,7,8},{9,10,11,12},{13,14,15,16}}); System.out.println(idenTest.dot(I,addr)); System.out.println(I.dot(idenTest,addr)); Addr<Double> addr1 = new Addr<Double>() { @Override public Double zero() { return 0.0; } @Override public Double one() { return 1.0; } @Override public Double add(Double lhs, Double rhs) { return lhs + rhs; } @Override public Double sub(Double lhs, Double rhs) { return lhs - rhs; } @Override public Double mul(Double k, Double hs) { return k * hs; } @Override public Double square(Double hs) { return hs * hs; } @Override public double prescription(Double hs) { return Math.sqrt(hs); } @Override public double div(Double hs, double nhs) { return hs / nhs; } @Override public Double div(Double hs, Double nhs) { return hs / nhs; } @Override public double acos(double hs) { return Math.acos(hs); } @Override public double toDegree(double hs) { return Math.toDegrees(hs); } @Override public int compair(Double lhs, Double rhs) { if (lhs - rhs > 0) { return 1; }else if (lhs - rhs < 0) { return -1; }else { return 0; } } @Override public boolean equals(Double lhs, Double rhs) { double dis = 1e-8; if (Math.abs(lhs - rhs) < dis) { return true; } return false; } @Override public Double neg(Double hs) { return -hs; } }; Ranks<Double> A = new Matrix<>(new Double[][]{{1.0,2.0},{3.0,4.0}}); System.out.println(A.invFunc(addr1)); System.out.println(A.invMatrix(addr1)); } }
运行结果
Matrix{values=[[1, 2, 3],[3, 4, 5]]}
[2, 3]
2
3
6
4
1
Vector{values=[3, 4, 5]}
Vector{values=[3, 5]}
Matrix{values=[[1, 3],[2, 4],[3, 5]]}
Vector{values=[366, 644]}
Matrix{values=[[155, 366],[275, 644]]}
Matrix{values=[[8, 10, 12],[13, 24, 35]]}
Matrix{values=[[6, 6, 6],[7, 16, 25]]}
Matrix{values=[[0.5, 1.0, 1.5],[1.5, 2.0, 2.5]]}
Matrix{values=[[1, 2, 3],[3, 4, 5]]}
Matrix{values=[[-1, -2, -3],[-3, -4, -5]]}
Matrix{values=[[3, 6, 9],[9, 12, 15]]}
Matrix{values=[[0, 0, 0],[0, 0, 0]]}
Matrix{values=[[1, 0, 0, 0],[0, 1, 0, 0],[0, 0, 1, 0],[0, 0, 0, 1]]}
Matrix{values=[[1, 2, 3, 4],[5, 6, 7, 8],[9, 10, 11, 12],[13, 14, 15, 16]]}
Matrix{values=[[1, 2, 3, 4],[5, 6, 7, 8],[9, 10, 11, 12],[13, 14, 15, 16]]}
Matrix{values=[[-1.9999999999999996, 0.9999999999999998],[1.4999999999999998, -0.49999999999999994]]}
Matrix{values=[[-1.9999999999999996, 0.9999999999999998],[1.4999999999999998, -0.4999999999999999]]}
Python代码
将原线性系统类修改为LinearSystemFunc
from .Vector import Vector from .Matrix import Matrix from ._global import is_zero class LinearSystemFunc: def __init__(self, A, b): assert A.row_num() == len(b), "矩阵的行数必须等于向量的长度" self._matrix_row = A.row_num() self._matrix_col = A.col_num() if isinstance(b, Vector): self.Ab = [Vector(A.row_vector(i).underlying_list() + [b[i]]) for i in range(self._matrix_row)] if isinstance(b, Matrix): self.Ab = [Vector(A.row_vector(i).underlying_list() + b.row_vector(i).underlying_list()) for i in range(self._matrix_row)] self.pivots = [] def _max_row(self, index_i, index_j, n): # 主元最大值所在的行 best, ret = self.Ab[index_i][index_j], index_i for i in range(index_i + 1, n): if self.Ab[i][index_j] > best: best, ret = self.Ab[i][index_j], i return ret def _forward(self): # 前向 i, k = 0, 0 while i < self._matrix_row and k < self._matrix_col: # 看Ab[i][k]位置是否可以为主元 max_row = self._max_row(i, k, self._matrix_row) self.Ab[i], self.Ab[max_row] = self.Ab[max_row], self.Ab[i] if is_zero(self.Ab[i][k]): k += 1 else: # 将主元归为一 self.Ab[i] = self.Ab[i] / self.Ab[i][k] for j in range(i + 1, self._matrix_row): self.Ab[j] = self.Ab[j] - self.Ab[j][k] * self.Ab[i] self.pivots.append(k) i += 1 def _backward(self): # 后向 n = len(self.pivots) for i in range(n - 1, -1, -1): k = self.pivots[i] # Ab[i][k]为主元 for j in range(i - 1, -1, -1): self.Ab[j] = self.Ab[j] - self.Ab[j][k] * self.Ab[i] def gauss_jordan_elimination(self): # 高斯-约旦消元法，如果有解，返回True;如果没有解，返回False self._forward() self._backward() for i in range(len(self.pivots), self._matrix_row): if not is_zero(self.Ab[i][-1]): return False return True def fancy_print(self): # 打印增广矩阵 for i in range(self._matrix_row): print(" ".join(str(self.Ab[i][j]) for j in range(self._matrix_col)), end=" ") print("|", self.Ab[i][-1])
新增加一个初等矩阵变换线性系统LinearSystemMatrix
from .Vector import Vector from .Matrix import Matrix from ._global import is_zero class LinearSystemMatrix: def __init__(self, A, b): assert A.row_num() == len(b), "矩阵的行数必须等于向量的长度" self._matrix_row = A.row_num() self._matrix_col = A.col_num() if isinstance(b, Vector): self.Ab = Matrix([Vector(A.row_vector(i).underlying_list() + [b[i]]) for i in range(self._matrix_row)]) if isinstance(b, Matrix): self.Ab = Matrix([Vector(A.row_vector(i).underlying_list() + b.row_vector(i).underlying_list()) for i in range(self._matrix_row)]) self.pivots = [] def _max_row(self, index_i, index_j, n): # 主元最大值所在的行 best, ret = self.Ab[(index_i, index_j)], index_i for i in range(index_i + 1, n): if self.Ab[(i, index_j)] > best: best, ret = self.Ab[(i, index_j)], i return ret def _forward(self): # 前向 i, k = 0, 0 while i < self._matrix_row and k < self._matrix_col: # 看Ab[i][k]位置是否可以为主元 max_row = self._max_row(i, k, self._matrix_row) self.Ab = self._change_matrix(i, max_row).dot(self.Ab) if is_zero(self.Ab[(i, k)]): k += 1 else: # 将主元归为一 self.Ab = self._mul_matrix(i, 1 / self.Ab[(i, k)]).dot(self.Ab) for j in range(i + 1, self._matrix_row): self.Ab = self._add_p_matrix(j, i, -self.Ab[(j, k)]).dot(self.Ab) self.pivots.append(k) i += 1 def _backward(self): # 后向 n = len(self.pivots) for i in range(n - 1, -1, -1): k = self.pivots[i] # Ab[i][k]为主元 for j in range(i - 1, -1, -1): self.Ab = self._add_p_matrix(j, i, -self.Ab[(j, k)]).dot(self.Ab) def _change_matrix(self, i, j): # 交换初等矩阵 identity = Matrix.identity(self._matrix_row) identity[(i, i)] = 0 identity[(j, i)] = 1 identity[(j, j)] = 0 identity[(i, j)] = 1 return identity def _mul_matrix(self, i, value): # 乘除初等矩阵 identity = Matrix.identity(self._matrix_row) identity[(i, i)] = value return identity def _add_p_matrix(self, i, j, value): # 加减初等矩阵 identity = Matrix.identity(self._matrix_row) identity[(i, j)] = value return identity def gauss_jordan_elimination(self): # 高斯-约旦消元法，如果有解，返回True;如果没有解，返回False self._forward() self._backward() for i in range(len(self.pivots), self._matrix_row): if not is_zero(self.Ab[(i, -1)]): return False return True def fancy_print(self): # 打印增广矩阵 for i in range(self._matrix_row): print(" ".join(str(self.Ab[(i, j)]) for j in range(self._matrix_col)), end=" ") print("|", self.Ab[(i, -1)])
矩阵类的inv方法也变成两个
from .Vector import Vector class Matrix: def __init__(self, list2d): self._values = [row[:] for row in list2d] @classmethod def zero(cls, r, c): # 返回一个r行c列的零矩阵 return cls([[0] * c for _ in range(r)]) @classmethod def identity(cls, n): # 返回一个n行n列的单位矩阵 m = [[0] * n for _ in range(n)] for i in range(n): m[i][i] = 1 return cls(m) def trans(self): # 返回矩阵的转置矩阵 return Matrix([[e for e in self.col_vector(i)] for i in range(self.col_num())]) def inv_func(self): # 返回矩阵的逆 from .LinearSystemFunc import LinearSystemFunc if self.row_num() != self.col_num(): return None n = self.row_num() ls = LinearSystemFunc(self, Matrix.identity(n)) if not ls.gauss_jordan_elimination(): return None invA = [[row[i] for i in range(n, 2 * n)] for row in ls.Ab] return Matrix(invA) def inv_matrix(self): # 返回矩阵的逆 from .LinearSystemMatrix import LinearSystemMatrix if self.row_num() != self.col_num(): return None n = self.row_num() ls = LinearSystemMatrix(self, Matrix.identity(n)) if not ls.gauss_jordan_elimination(): return None invA = [[row[i] for i in range(n, 2 * n)] for row