线性代数计算的性能
压缩数组
建立 Wolfram 语言的挑战之一是确保系统具有足够的一般性来支持符号计算并具有足够快的速度来进行有效的数值计算. 关于这些及其它目标在 "Wolfram 语言的设计原理" 一节中介绍. Wolfram 语言提供的一项技术是对 Wolfram 语言表达式的重要类型使用专业化的表示. 对于线性代数,这涉及到版本 4.0 中首次提出的所谓压缩数组技术.
vec = Table[Random[], {10000}];
vec//ByteCountDo[vec ^ 2, {50}]//Timingvec[[1]] = xvec//ByteCountDo[vec ^ 2, {50}]//Timing第一个向量使用的内存较少,支持更快的计算,原因是它是由压缩数组构建的. 压缩数组代表的是机器整数、机器实数以及整数和虚数部分均为机器实数的复数的矩形张量.
vec = Table[Random[], {1000}];
Developer`PackedArrayQ[ vec]vec[[2]] = 10.1;
Developer`PackedArrayQ[vec]vec1 = vec ^ 2;
Developer`PackedArrayQ[vec1]vec[[2]] = x;
Developer`PackedArrayQ[vec]vec = Table[Random[], {1000}];
vec[[2]] = 0;
Developer`PackedArrayQ[ vec]压缩数组函数
| Developer`PackedArrayQ[expr] | 若 expr 是压缩数组则返回 True |
| Developer`ToPackedArray[expr] | 尽可能将 expr 转化为压缩数组 |
| Developer`ToPackedArray[expr,t] | 尽可能将 expr 转化为类型为 t 的压缩数组 |
| Developer`FromPackedArray[expr] | 将压缩数组转化为 expr |
| Developer`PackedArrayForm[expr] | 以摘要形式打印 expr 的压缩数组 |
pack = Developer`ToPackedArray[{1, 2, 3}];
Developer`PackedArrayQ[ pack]pack = Developer`ToPackedArray[{1, 2, 3.5}];
Developer`PackedArrayQ[ pack]pack = Developer`ToPackedArray[{1, 2, 3.5}, Real];
Developer`PackedArrayQ[ pack]Developer`PackedArrayForm[pack]pack = Developer`ToPackedArray[{{1, 2}, 2}];
Developer`PackedArrayQ[ pack]压缩数组运算
pack = Developer`ToPackedArray[{1, 2, 3}];
pack === {1, 2, 3}MatchQ[{1, 2, 3}, _List]MatchQ[pack, _List]与非压缩列表相比,许多函数在压缩数组的应用更为有效,但结果将相同,不管所用的是压缩或非压缩形式. 为了确保前后一致,有时系统必须要将压缩数组进行解压缩. 由于这将是导致低效的潜在原因,可有一则信息来告知用户什么已被解压缩. 通过一个系统选项启用这项功能.
SetSystemOptions[ "PackedArrayOptions" -> {"UnpackMessage" -> True}];
v = Range[100];
Developer`PackedArrayQ[ v]v[[1]] = 5.6;SetSystemOptions[ "PackedArrayOptions" -> {"UnpackMessage" -> False}];压缩数组总结
mat = {{1., 2.}, {4., 5.}};
Developer`PackedArrayQ[ mat.mat]通常很难利用诸如 Developer`ToPackedArray 等函数来启用压缩数组. 主要问题是确保矩阵中各元素的一致性,例如确保不要将整数和实数混合.
编程效率
本节讨论如何编写高效的 Wolfram 语言程序来求解线性代数问题.
性能的度量
如果您想研究自己编写的 Wolfram 语言程序代码的性能,有一些与时间相关的函数可以使用.
| Timing[expr] | 计算 expr 并返回 CPU 所需时间和结果的列表 |
| AbsoluteTiming[expr] | 给出所用的绝对时间,计算 expr |
mat = Table[Random[], {1000}, {1000}];
Timing[mat.mat;]另一个有用的函数是 TimeConstrained.
| TimeConstrained[expr,t] | 尽量计算 expr,在 t 秒后停止计算 |
| TimeConstrained[expr,t,failexpr] | 如果没有达到时间限制,返回 failexpr |
mat = Table[Random[], {1000}, {1000}];
TimeConstrained[ Do[mat.mat, {10}], 1]循环向量化
常见的一种运算涉及到在矩阵的元素上进行某种类型的迭代. 由于矩阵可能很大,效率显然是至关重要的.
编写高效代码最重要的方式之一就是要避免显式部分引用,特别是在内部循环中. 这是用户代码效率低下的主要来源,使用同时作用在矩阵或向量的所有元素上的函数可以避免这种问题. 这些向量化的操作允许 Wolfram 语言使用高度优化的代码,运行速度更快.
创建列表
这个问题涉及到创建一个列表,然后用一定值进行初始化. 例如,计算
v〚i〛 ⟹ Sin[i]
一种方法是创建列表然后使用一个 For 循环在各个元素上进行迭代从而对这些元素进行赋值.
n = 10;
v = Table[0, {n}];
For[i = 1, i <= n, i++, v[[i]] = Sin[2.0 Pi i]];
n = 10;
v = Table[ Sin[2.0 Pi i], {i, 1, n}]
一种更为快捷的方式是使用优化函数 Range 来创建列表,然后使用向量化的运算.
n = 10;
v = Sin[2.0 Pi Range[1., n]]
用于创建列表和矩阵的诸如 Table 和 Range 等高度优化的命令将在 "构建矩阵" 中讨论.
更新列表
这个问题涉及到列表元素的更新. 例如,计算每个元素的 Sin.
v〚i〛 ⟹ Sin[v〚i〛]
一种方法涉及在列表的每个元素上进行迭代,将每个元素用更新后的元素替换.
Do[v[[i]] = Sin[v[[i]]], {i, n}];
v = Sin[v]
v〚i〛 ⟹ v〚i〛/i^2
解决这个问题的一种方法涉及到在列表的元素上进行迭代,将每个元素除以索引的平方.
n = Length[v];
Do[ v[[i]] = v[[i]]/i^2 , {i, n}];
d=Range[n]^2;
v=v/d
可列表运算的使用在"可列表性" 中讨论.
使用内置支持
Wolfram 语言具有范围非常广泛的函数,包括许多在矩阵计算中有重要应用的列表处理函数. 在开始工作之前,查看一下是否有内置函数能够帮您实现目的是很有必要的. 例如,如果所进行的工作涉及卷积或相关,那么应该使用适当的函数,这些函数几乎肯定会更快.
另外要注意的重要一点是,Wolfram 语言能够表示范围广泛的数学对象. 例如,Wolfram 语言能用于解决实际的线性方程问题,这些线性方程由矩阵来表示. 有时,方程形式的表述性更强也更方便,从而导致效率提高. 这个主题将在"将方程转化为稀疏矩阵"中讨论.
列表的结构化操纵是应用于许多线性代数计算的内置函数的一个领域. 诸如 Part, Take, Drop, Transpose, PadLeft, PadRight, RotateLeft 和 RotateRight 等 Wolfram 语言函数是经过高度优化的,并有众多应用. 有关于此的更多细节将在"结构化运算"中讨论. 如果您可以编写自己的计算以使用这些函数,通常将是非常理想的.
现在给出一个例子来讨论通过平铺复制来拓展矩阵的不同方法.例如,通过重复填充输入矩阵
一种缓慢的方法
extendMatSlow[mat_, nRow_, nCol_] :=
Module[ {m, n, mf, nf, final},
{m, n} = Dimensions[mat];
mf = m * nRow;
nf = n * nCol;
final = Table[0, {mf}, {nf}];
Do[
Do[final[[m * (i - 1) + i1, n * (j - 1) + j1]] = mat[[i1, j1]], {i1, m}, {j1, n}], {i, nRow}, {j, nCol}];
final
]mat = {{1, 2, 3}, {4, 5, 6}};
extendMatSlow[mat, 2, 3]//MatrixFormmat = Table[Random[], {100}, {20}];
extendMatSlow[mat, 10, 15];//Timing这种方法存在几种速度上的缺陷,最严重的是在部分索引上进行内循环. 如果替代以向量化代码,速度将大大加快. 在"性能的度量"中已经解释过在内部循环中进行显式部分引用,将导致代码的低效.
sp = SparseArray[{{1, 2} -> 4, {2, 3} -> 1}, {2, 3}];
extendMatSlow[sp, 2, 3]一种较快的方法
extendMatFaster[mat_, nRow_, nCol_] :=
Module[ {m, n, mind, nind},
{m, n} = Dimensions[mat];
mind = Range[m];
mind = Flatten[Table[mind, {nRow}]];
nind = Range[n];
nind = Flatten[Table[nind, {nCol}]];
mat[[mind, nind]]
]mat = {{1, 2, 3}, {4, 5, 6}};
extendMatFaster[mat, 2, 4]mat = Table[Random, {100}, {20}];
extendMatFaster[mat, 10, 15];//Timing这种方法要快得多;它使用函数的 Part 非常灵活并经过了高度优化. 这将在"得到矩阵的一部分"中进一步介绍.
sp = SparseArray[{{1, 2} -> 4, {2, 3} -> 1}, {2, 3}];
extendMatFaster[sp, 2, 3]速度快并且更简洁
mat = {{1, 2, 3}, {4, 5, 6}};
PadRight[mat, {4, 9}, mat]mat = Table[Random, {100}, {20}];
PadRight[mat, {1000, 300}, mat];//Timing这种方法也更加简洁,并且对于小型输入矩阵来说可能是最快的. 对于大型矩阵,其速度与已经介绍过的 Part 范例中的速度相似. PadRight 的用法在"矩阵的扩展"中介绍.
sp = SparseArray[{{1, 2} -> 4, {2, 3} -> 1}, {2, 3}];
PadRight[sp, {4, 9}, sp]矩阵内容
Wolfram 语言中的矩阵可以由 Wolfram 语言所持有的所有不同类型对象构建. 这些矩阵可以包含机器精度浮点数、任意精度浮点数、复数浮点数、整数、有理数以及一般的符号量. 这些不同类型在"矩阵类型"中讨论.
Wolfram 语言可以操作这些不同类型的矩阵,这一点给系统强大的力量. 但是,它对数值浮点计算的效率有一些影响:首先来自于符号和数值项的组合;其次来自于不同类型数字的组合;第三来自于整数矩阵.
混合的符号/数值矩阵
如果在一个矩阵中将符号型和数值型元素混合,线性代数函数将如"混合模式矩阵"中所讨论的,把该矩阵视为符号矩阵. 应该注意的是,操作符号矩阵的方法通常比操作数值矩阵要慢,因为对于数值矩阵,有许多高度优化的技术可用.
mat = Table[Random[], {200}, {200}];
mat.mat;//Timingmat[[2, 2]] = x;
mat.mat;//Timing混合数值类型矩阵
如果矩阵中混合了不同的数值类型,线性代数函数将把矩阵视为数值型,并将矩阵强制转换为最低的精确度. 这在"混合模式矩阵"中讨论. 对于线性代数函数,这样做对性能所带来的影响将不会像处理符号矩阵那么严重(这一点在前面一节中已经说明).
mat = Table[Random[], {200}, {200}];
mat.mat;//Timingmat[[2, 2]] = 1;
mat.mat;//Timingn = 300;
vec = Table[0, {n}];
Do[ind = Random[Integer, {1, n}];vec[[ind]] = Total[1.5 + 2 vec], {100}];//TimingDeveloper`PackedArrayQ[vec]n = 300;
vec = Table[0.0, {n}];
Do[ind = Random[Integer, {1, n}];vec[[ind]] = Total[1.5 + 2 vec], {100}];//TimingDeveloper`PackedArrayQ[vec]这个问题的解决办法是,如果将要进行实数的操作,一定要确保使用实数进行初始化.
整数矩阵
mat = {{1, 2}, {3, 4}};
Do[ Eigenvalues[mat], {100}];//Timingmat = {{1, 2}, {3, 4}};
mat = N[mat];
Do[ Eigenvalues[mat], {100}];//TimingEigenvalues[{{1, 2}, {3, 4}}]Eigenvalues[N[{{1, 2}, {3, 4}}]]当然,Wolfram 语言的一个强项是,它能够对涉及整数矩阵的计算返回精确结果. 但应该注意的是,许多纯数值系统将返回近似结果.
这类问题的解决方法是,如果进行的是关于近似实数矩阵的操作,应该以近似实数开始.
表达式效率
"用作 Wolfram 语言表达式的矩阵" 对 Wolfram 语言用表达式——贯穿整个系统的一种统一的数据类型——的形式表示矩阵做出了解释. 它介绍了系统使用这种统一格式的优越性. 本节在此讨论的是使用 Wolfram 语言表达式对效率带来的一些影响.
矩阵的更新
vec1 = {1, 2, 3};
vec2 = vec1;vec1[[1]] = 4;
vec1vec2vec = Range[1., 50000];
Do[vec[[i]] = i - 1, {i, 1., 50000}];//Timingvec = Range[1., 50000];
Do[vec1 = vec;
vec[[i]] = i - 1, {i, 1., 50000}];//Timingvec = Range[1., 50000];
Do[ReplacePart[vec, N[i] - 1, i], {i, 1, 50000}];//Timing使用的是先前介绍的向量化运算,能够避免在矩阵上进行迭代,并将对运算的复制次数最小化. 如果在矩阵上的迭代无法避免,一个明智的方法是尽可能的使循环简洁从而避免额外复制矩阵.
整体过程所需时间与在每一步采用复制方法的列表长度的平方成比例. 如果不进行复制,则时间分配是线性的. 输入列表越长,影响越大,并且代码的性能将开始降低.
矩阵的增补
在 Wolfram 语言中,表达式按数组方式执行. 这允许快速访问表达式的任何具体内容. 由于矩阵是正规的 Wolfram 语言 表达式,它们也是数组. 虽然数组元素的访问速度非常快,添加元素却较慢. 这通常需要一个副本,这种情况应该避免. 后面的例子说明了增加元素的代价.
vec = {};
Do[AppendTo[vec, i], {i, 10000}];//Timingvec = Table[i, {i, 10000}];//Timingnum = 10000;
list = Table[Random[], {num}];res1 = {};
Map[If[# > 0.5, AppendTo[res1, #]]&, list];//Timingres2 = Flatten[Map[If[# > 0.5, #, {}]&, list]];//Timingres3 = Reap[Map[If[# > 0.5, Sow[#]]&, list]][[2, 1]];//Timingres1 == res2 == res3整体过程所需时间与在每一步采用复制方法的列表长度的平方成比例. 如果不进行复制,则时间分配是线性的. 输入列表越长,影响越大,并且代码的性能将开始降低.
相关技术笔记
-
▪
- Wolfram 语言中的线性代数 ▪
- 矩阵与张量运算 ▪
- 稀疏数组的运用 ▪
- 矩阵计算 ▪
- 矩阵类型
