大数据计算方法 课件 第1-4章 绪论、数据分析基础-稀疏矩阵解法

1大数据计算方法Chap1：绪论主要内容1.绪论

(误差分析/Matlab基础/浮点算术体系)2.数据分析方法概述(数据分析的基本概念与模型的评价方法)3.线性方程组求解与矩阵分解(LU/Cholesky分解/QR分解与最小二乘)4.大型稀疏线性方程组求解(稀疏矩阵/直接解法/CG/GMRES算法)5.矩阵的特征值与奇异值(QR算法/Krylov子空间迭代法/SVD算法)6.优化问题与凸优化

(优化问题基础/凸优化/举例:几何优化)7.数据分析中的凸优化(最大间隔分类/正则化回归/压缩感知)8.非线性方程与无约束优化(牛顿法/简单方法/基于求导的方法)9.约束优化问题解法与对偶原理

(等式约束/不等式约束/对偶原理)10.基于矩阵分解的数据挖掘(PCA及应用/随机化算法/非负矩阵分解)大数据计算方法2主要借助Matlab学习(也涉及Python)Matlab的替代引言Matlab与数值计算软件浮点算术体系与误差分析Outline3引言概览“十大算法”问题、以及求解步骤误差分析基本概念回顾41950年之前算法线性插值牛顿法高斯消去法工具手算数学手册数据表格数值计算的历史51950年代之后现代计算机的发明催生了大量

新的数值计算方法/算法数值计算的历史6数值计算：主要研究求解连续数学问题的算法（其中包含、但不局限于计算误差的研究）Evolutionofnumericalcomputingfromothersciencesandengineeringdisciplines数值计算是求解很多实际工程问题的重要工具数值计算8PageRank算法对1010x1010、甚至更大的矩阵计算其特征向量——世界上最大规模的矩阵计算航空公司使用优化算法来确定机票

价格、航班调度和航油供应更多现实世界的例子数值计算9汽车公司通过对汽车碰撞的计算机模拟来发现潜在的安全隐患金融公司通过统计仿真工具来预测股票价格常微分方程(Ordinarydifferentialequation,ODE)电路的瞬态分析(transientanalysis)数值计算问题举例10

...代数方程偏微分方程(Partialdifferentialequation,PDE)热分析(thermalanalysis)：求温度分布数值计算问题举例11比重比热容热导率热源的功率密度热流(单位面积热通量)温度发热功率密度(傅里叶热传导定律)线性回归(Linearregression)也称为响应面建模(responsesurfacemodeling)数值计算问题举例12求解超定线性方程组来得到参数a和b(方程数目比变量多)带正则项/约束的线性回归(regularization)数值计算问题举例13最小二乘误差L1-范数正则项原始图片上采样

恢复出的图片分类问题(Classification)最大间隔分类-支撑向量机(SVM)方法数值计算问题举例14求解凸二次规划得到W和CL2-范数(最大间隔)WenjianYu1515JackDongarra获2021年图灵奖主要源于其在数值算法和库方面的开创性贡献，正是他的研究使高性能计算软件在四十多年里与硬件的指数级改进保持同步。Toptenalgorithmsofthecentury1.1946LosAlamos国家实验室的J.vonNeumann,S.Ulam和N.Metropolis发展出来的Metropolis算法

（属于MonteCarlo方法；拒绝采样/MCMC,生成给定概率分布随机点)2.1947RAND公司的G.Dantzig创造的解线性规划的单纯型算法(simplexmethod)3.1950UCLA与美国国家标准局数值分析所的M.Hestenes,E.Stiefel和C.

Lanczos开创的Krylov子空间迭代法(CG算法，Lanczos过程)4.1950’s矩阵分解方法,由OakRidge国家实验室的A.Householder引入数值线性代数中(矩阵计算研究中掀起革命)“Wetriedtoassemblethe10algorithmswiththegreatestinfluenceonthedevelopmentandpracticeofscienceandengineeringinthe20thcentury”Editors16通用性,统一,reuse,舍入误差分析,updating,软件开发Toptenalgorithmsofthecentury5.1957IBM的J.

Backus领导的小组开发的Fortran优化编译器6.1959-61

伦敦FerrantiLtd.的J.G.F.Francis发明的QR算法，稳定地计算中、小规模矩阵的所有特征值7.1962伦敦ElliottBrothers,

Ltd.的TonyHoare提出快速排序算法(Quicksort)8.1965IBMWatson研究中心的J.Cooley与Princeton大学及AT&TBell实验室的J.Tukey共同提出了的FFT算法9.1977BrighamYoung大学的H.

Ferguson和R.Forcade提出的整数关系侦察算法(实验数学,简化量子场理论的计算)10.1987Yale大学的L.Greengard和V.

Rokhlin发明了快速多极算法(fastmultipolealgorithm,多体问题,线性复杂度)大多数都属于或涉及数值计算的范畴!17Newton/quasi-Newton法(QZ/关于Schur,SVD)PagerankJPEGKalmanfilterN.Higham,ThePrincetonCompaniontoAppliedMathematics,PrincetonUniversityPress,2015数值模拟的步骤建立数学模型（需要相关学科背景）研究数值求解方程或最优化的算法通过计算机程序(软件)实现算法在计算机上运行程序(软件),进行数值实验数值计算的问题与步骤数值计算（科学计算、计算模拟、数值模拟）没有解析解的问题：虽有解析表达式(解)，但无法(难以)计算:模拟通常难以达到的实验条件(时间、金钱成本)天体物理研究汽车碰撞实验，芯片流片(制造)前的性能仿真除了解方程的问题，还有更多的优化问题18sin(x)

数值计算的问题数值模拟的步骤（续）将计算结果用较直观的方式输出，如图形可视化方法解释和验证计算结果。若需要，重复上面的某些步骤上述各步骤相互间紧密地关联，影响着最终的计算结果和效率（问题的实际应用背景也左右着方法的选择）重视计算中的扰动原始数据、参数的扰动计算过程中的扰动扰动可能遵循一定的规律19它们对结果如何影响？(UncertaintyQuantification)(Round-offErrorAnalysis)近似造成误差模型误差：建模时的近似数据误差：测量等原因造成误差、前一步计算的结果截断误差(方法误差)：数学上的近似舍入误差：实数用有限位来表示(数值计算中总存在)误差计算地球的表面积：

20模型误差数据误差数据误差舍入误差截断误差与舍入误差统称“计算误差”

误差

数据误差引起的

f(x)=sinx,x=121

步长h

(截断误差、舍入误差)大体上

问题敏感性与算法稳定性

病态问题?22(相对条件数)Matlab与数值软件准备知识黄金分割比斐波那契数分形蕨幻方数值软件23MATLAB简介(1)MATrixLABoratory(矩阵实验室),

集成可视化快速原型编程工具,数值/矩阵计算功能强(先进算法)大量专题工具箱(Toolbox),为专业应用提供便利Matlab(作为编程语言)C,C++,Fortran第四代编程语言第三代编程语言

编译方式解释器，或JIT加速器(v6.5后)编译器

申明变量？不需要需要

开发时间较快较慢

运行时间较慢较快

开发环境集成环境(编辑器、调试器、命令历史、变量空间、profiler、编译器)--24MATLAB2020b的函数手册>>helpsvdsketchMATLAB简介(2)命令窗口，将Matlab作为计算器使用变量空间命令历史其他窗口25MATLAB简介(3)程序编辑窗口注释语句%帮助窗口26找不到Matlab软件怎么办？尝试GNUOctave

/0.6180…，可能是世界上最有趣的数？Matlab功能表示、显示数的格式数值、符号(解析)求解多项式方程符号变量转为数值函数句柄与函数的绘制脚本程序与绘图功能一个较复杂的函数程序(截断n项的连分数)黄金分割比(1)27

黄金分割比(2)>>phi=(1+sqrt(5))/2>>formatlong

28求解多项式方程用降次排列的系数表示多项式，再用roots命令数值求解用SymbolicMathToolbox解析求解命令solve命解方程:解析(符号)求解运算长无法求解复杂方程符号变量转化为数值vpa,double,single

命令注意：抑制输出符“;”

黄金分割比(3)>>p=[1-1-1];>>r=roots(p)

>>vpa(phi,50)>>phi=double(phi)29>>symsx;r=solve(x-1==1/x)>>phi=r(2)

黄金分割比(4)

>>f=@(x)x-1-1./x>>phi=fzero(f,1)>>holdon>>plot(phi,0,'o')

30>>fplot(f,[0,4])脚本程序与绘图功能Matlab程序.m文件有两种:脚本(script),函数(function)一个脚本程序实现绘制黄金矩形的功能程序内容为：黄金分割比(5)>>editgoldrect.m%GOLDRECTGoldenRectangle…plot(x,y,'b',u,v,'b--')text(phi/2,1.05,'\phi')…axisequalaxisoffset(gcf,'color','white')注释、plot,text,axis,set命令getcurrentfigure%遵循真实尺寸%删除坐标轴31

黄金分割比(6)(类似于C语言中的system)

32斐波那契序列(LeonardoP.Fibonacci与他的兔子问题)f1=1,f2=2,fn=fn-1+fn-2第一个函数fibseq.m，计算出前n个斐波那契数生成全零向量函数:zeros第二个函数fibn.m，输出第n个斐波那契数recursive(递归):

递归函数形式简洁、但运行开销大!计时命令：tic,toc相邻斐波那契数之比~截断的连分数~黄金分割比斐波那契数>>f=zeros(n,1)>>tic,fibn(35),toc>>tic,f=fibseq(35);f(35),toc>>f=fibseq(40);>>f(2:40)./f(1:39)33试试fibn(60)?逐项运算符./“神奇”的通式原理平面上精心策划的点集合、无限的随机生成点的过程

,有4种变换,第1种是每次随机地应用4种变换中的一种，生成新的点无限生长程序:fern含n个点的函数:finitefernMatlab命令程序中含较复杂的绘图命令定义矩阵:‘;’while语句rand命令生成0~1之间均匀随机数用imread,image可看图片分形蕨>>fern>>finitefern(100000)>>A1=[.85.04;-.04.85];>>r=rand

>>F=imread('fern.png');>>image(F)34特殊的矩阵~magic命令《周易》中“洛书”(3阶幻方)的传说行,列,正/反对角线元素之和均相等Matlab命令验证magic(3)是幻方这些和等于多少?三阶幻方矩阵有哪几种?幻方(1)>>A=magic(3)>>sum(A)>>sum(A')'>>sum(diag(A))>>sum(diag(flipud(A)))fork=0:3rot90(A,k)rot90(A',k)end35>>det(A)>>formatrat>>X=inv(A)%逆时针转矩阵的行列式/逆矩阵三阶幻方看看矩阵的范数、

特征值更高阶幻方AlbrechtDurer的版画MelancoliaII4阶幻方有880个

(5阶~275305224个!)矩阵的秩生成幻方的算法幻方矩阵的三维表面图幻方(2)>>loaddurer>>image(X)>>colormap(map)>>axisimage>>inv(A)>>

rank(A)(丢勒)36>>norm(A)>>eig(A)>>editmagic.msurf(magic(n));axisoff;n阶幻方的最大特征值是？(.mat文件)常数:pi,i,j算术运算:+,-,*,/,^,.*,./,.^逻辑关系:==,~=,>,>=,<,<=逻辑运算:&,|,~,&&,||帮助:help,doc,lookfor,which,demo输出:disp,sprintf,;,format初等函数:sin,cos,tan,sinh,asin,exp,log,log10,sqrt

变量记录:whos,clear,save,load,ans,diary向量:[...,...],[...;...],length矩阵:[...,...;...,...],ones,zeros,eye,rand,size,diag,tril,triu绘图:plot,subplot,loglog,fplot,hold,plot3,figure,close文件:edit,type,ls,path编程:function,if,for,while,end,@退出:exit一些常用的matlab命令37数值软件/程序包数值计算的软件与程序包解决一些共同问题，促进各个科研领域的工作了解基本原理，学习算法设计和实现技巧成为聪明的软件/程序包使用者存在形式和资源免费(netlib,github,TOMS…)付费(NR,…)Fortran,C,C++,Matlab,Python,Julia,…源代码使用，或API调用;交互式集成环境的软件38数值软件/程序包矩阵计算有关程序包LAPACK,/lapack/index.htmlBasiclinearalgebrasubprogram(BLAS)

/blas/厂家提供的BLAS;ATLAS,GotoBLAS,OpenBLAS*ScaLAPACK:Distributed-memoryvariantMAGMA:MatrixAlgebraforGPU&MulticoreArch.更高层的程序库与环境Intel公司MathKernelLibrary(inC)Eigen(aC++templatelibrary):Matlab,Octave,Python(Numpy),…392016DenseLinearAlgebraSoftwarePackagesSurvey

(J.Dongarra)87%40%22%No.1No.2(有20%需要自己写线性代数程序)/lapack/lawnspdf/lawn290.pdf“asmuchaspossiblecomputationisperformedbycallstoBLAS”“exploitLevel3BLAS”Linux:95%,Mac:32%,Win:22%(LAPACKE)浮点算术体系与误差计算机中的浮点数浮点数与舍入误差数值计算中如何减小误差40计算机中的浮点数

>>realminans=2.2251e-308>>realmaxans=1.7977e+308realmin=2-1022realmax=(2-eps)2102341[-1022,1023]

缺哪些数?(非规范数)计算机中的浮点数

浮点数系统IEEE单精度23-1261275.960

10-8IEEE双精度52-102210231.110

10-16

例:一个简单浮点数系统,

42p=2,L=-1,U=1浮点数与舍入误差

,,>>formathex;t>>a=4/3>>b=a-1>>c=3*b>>e=1-c=eps/2

1.110

10-16

(double型)43E+1023尾数f舍入——“四舍五入”

Matlab用浮点数表示整数浮点整数，只要”结果数的位数不太多”，其运算没有误差若加减乘除/平方根结果还是整数，无误差命令[F,E]=log2(x)实现浮点数的分拆:x=F

2E,F

[0.5,1)由于舍入误差，奇异矩阵在计算时可能并不奇异例如,inv(magic(4))

一个特殊的7次多项式Matlab中相关现象>>x=0.988:.0001:1.012;>>y=x.^7-7*x.^6+21*x.^5-35*x.^4+35*x.^3-21*x.^2+7*x-1;>>plot(x,y)(x-1)744尾数52位,

253的整数必可精确表示发生了严重的抵消现象

?~9.01015>>holdon;y2=(x-1).^7;plot(x,y2)

(反过来用pow2)抵消与“大数吃掉小数”

45

另一个例子一元二次方程求根公式的例子解为：计算x2时可能出现的抵消现象可类似地解决46

如何避免？数值计算中如何减小误差

总误差计算误差数据

传递误差截断误差舍入误差如何评估大小？分析具体的计算方法向后误差分析；区间分析法；一般很难定量分析问题敏感性(条件数);直接近似分析如何减小误差？选择截断误差小的算法选稳定的算法；减小舍入误差的建议；采用更高精度浮点数变换问题形式(计算过程)，改善敏感性4748大数据计算方法Chap2：数据分析方法概述基本概念分类回归分析聚类分析和降维Outline49数据分析基本概念性能度量实验测试50

基本概念51

基本概念52性能度量查准率（又称精确率，Precision）和查全率（又称召回率，Recall）查准率：衡量预测模型预测正确的正（负）样本个数占该模型所有分类为正（负）样本个数的比例查全率：模型预测正确（错误）的正样本个数占所有的正（负）样本个数的比例通过混淆矩阵进行描述和计算表2.1二分类的混淆矩阵基本概念53预测为“正”的样本数预测为“负”的样本数真实为“正”的样本数TruePositive(TP)FalseNegative(FN)真实为“负”的样本数FalsePositive(FP)TrueNegative(TN)

基本概念54图2.1P-R曲线示意图性能度量F1-score通过计算查全率和查准率的调和平均值来对两者进行综合评价实际数据集中，会出现样本类别分布不平衡的情况，P-R曲线会受到样本分布变化而产生较大的震荡，由此引入：ROC曲线（ReceiverOperatingCharacteristic）被试者在不同判断标准下所得的虚报概率FPR为横坐标，以击中概率TPR为纵坐标当数据集正负样本的分布发生变化或正负样本不平衡时，该曲线能够保持不变基本概念55

性能度量AreaUnderrocCurve(AUC)ROC曲线下的面积，度量分类模型好坏的量化指标AUC值越接近于1.0，该模型的性能越好不受样本分布变化和不平衡的影响等错误率(EqualErrorRate，EER)衡量分类模型好坏的客观评价指标ROC曲线上错误分类一个正样本或负样本概率相等时的点基本概念56图2.2

ROC曲线示意图实验测试模型的泛化性能需要在实验测试过程中进行评估数据集充足情况下，可以被随机划分为训练集、验证集和测试集训练集数据用于模型的学习，验证集数据用于模型选择，而测试集数据用于对最终模型泛化性能的评估数据的划分需要尽量满足分布的一致性以及互斥性常用的训练-测试数据集的划分方法包括：留出法、交叉验证法及自助法基本概念57实验测试留出法(Hold-out）直接按比例将数据集划分为三个互斥的子数据集常见的划分方案是将2：3或4：5比例的样本用于训练，其余样本则用于验证集和测试集的划分基本概念58图2.3

用留出法进行数据划分示意图实验测试交叉验证法K折交叉验证（K-foldCrossValidation）和留一交叉验证（Leave

One-Out）基本概念59

图2.4

K折交叉验证示意图K折交叉验证的一种特例，令K=样本数m

基本概念60分类二分类K近邻法线性判别分析逻辑回归61

分类62

分类63

分类64

K近邻法（K-NearestNeighbor,KNN）实质上是对训练样本特征空间的划分，其中K值的选择、距离度量方法及分类的决策规则是该算法的三个要素（1）K值：算法的一个超参数。当K值较小时，x的类别只与它非常接近的样本有关。学习的近似误差会减小，但估计误差会增大。如果K值较大，x的类别就由较大范围的样本决定，可以减少估计误差，但近似误差增加。分类65

分类66K近邻法（K-NearestNeighbor,KNN）（3）决策规则：往往采用多数表决的决策策略，即由待分类样本的K个最邻近的训练样本中的多数类决定其类别K近邻算法的复杂度取决于训练集样本的个数。当算法应用于样本较多、特征维度较高的大数据场景时非常耗时，必须采用高效的数据结构存储训练数据，以减少距离计算次数。KD（K-DimensionalTree）树是一种能够实现快速K近邻检索的方法。它是一种对K维空间中的样本点进行存储以便对其进行快速检索的树形数据结构，即二叉树。分类67K近邻法（KNN）KD树

实质：对

K

维空间进行划分构造方法与二叉树类似构造过程：不断用垂直于坐标轴的超平面来切分K维空间，得到一系列K维超矩形区域，而最终KD树的每一个节点就对应于一个K维超矩形区域，如图2.5所示。KD树搜索将搜索范围限制在一个局部超巨型体区域，不需要计算距离，大大减少了距离计算次数，算法的时间复杂度可降低至O(nlogn)。分类68图2.5KD树示意图

分类69线性判别分析线性判别分析（LinearDiscriminantAnalysis，LDA）是一种经典的线性分类模型主要思想：在训练阶段，算法找到一条直线，使得所有训练样本投影到该直线上后，可以让同类样本的投影点尽可能接近，而异类样本的投影点尽可能远离；在测试阶段，将待分类样本同样投影到该直线上，再根据投影点的位置来确定样本的类别,如图2.6所示。分类70图2.6LDA算法原理示意图

分类71

分类72

分类73

分类74

分类75

分类76

分类77

分类78回归分析线性回归分析非线性回归分析欠拟合和过拟合79回归分析是确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法按照涉及的因变量的多少，分为简单回归分析和多重回归分析按照自变量的多少，可分为一元回归分析和多元回归分析按照自变量和因变量之间的关系类型，可分为线性回归分析和非线性回归分析回归分析80

回归分析81

回归分析82

回归分析83

欠拟合和过拟合衡量模型的好坏不仅要求模型对训练数据集有很好的拟合（训练误差小），更希望它可以对未知数据集（测试集）有很好的拟合结果（泛化误差小）。而模型的过拟合（overfitting）和欠拟合（underfitting）都是导致模型的泛化能力不高的两种常见原因，其中过拟合是机器在学习中经常面临的棘手问题。不同复杂度的模型会得到不同的训练误差和泛化误差回归分析84欠拟合和过拟合下图中，以线性式拟合的模型用直线去逼近各个训练样本点，得到了较大的训练误差和测试误差，这种情况就是“欠拟合”。85解决欠拟合问题较容易，可以通过增加模型复杂度或特征维度等方式。欠拟合和过拟合八次多项式模型的逼近比三次多项式模型更接近训练样本点，但与实线表示的目标函数已经产生了背离，产生了较大的测试误差。十六次多项式模型得到了非常小的训练误差，但测试误差非常大。这两种模型都是在训练集上表现很好，而在测试集上表现很差，这种情况称为“过拟合”问题。86欠拟合和过拟合造成过拟合的原因：训练数据集样本单一，样本不足训练数据中噪声干扰过大。噪声指训练数据中的干扰数据。模型过于复杂（重要因素）。如模型太复杂，模型相当于可以做到“死记硬背”出训练数据的所有信息，但对没有见过的数据则不会变通，泛化能力太差。优化过于复杂模型过程就是寻找能够解释已知数据且最简单的模型的过程人们常采用正则化（Regularization）、早停（Early

stopping）、随机失活（Dropout）、数据增强等方法来抑制过拟合。回归分析87欠拟合和过拟合早停法：模型迭代训练过程中，在模型对训练数据集收敛之前就停止训练，以防止模型由于过度训练而产生过拟合。随机失活法：随机失活只应用于人工神经网络模型的过拟合抑制，它通过随机丢弃一些神经元连接来达到降低模型复杂度的目的。数据增强：让模型具有更好泛化能力的方法还可以从增加训练数据的角度。数据增强就是通过增加训练集的额外副本来增加训练集的大小，进而改进模型的泛化能力。回归分析88聚类分析和降维聚类算法的性能度量距离计算K均值聚类算法层次聚类密度聚类降维89聚类分析属于探索性的数据分析方法聚类(Clustering)

是针对大量数据,通过对无标记训练样本的学习来发现数据的内在性质及规律，再根据数据本身的特性将相似数据进行分组由于训练样本对应的类别是未知的，聚类是一种无监督学习方法聚类试图将数据集中的样本划分为若干个通常是不相交的子集，每个子集称为一个“簇”（cluster）聚类分析和降维90

聚类分析和降维91聚类算法的性能度量聚类应满足同一簇样本尽量彼此相近，而不同簇的样本则应尽可能不同分为两类指标：外部指标和内部指标外部指标是将聚类结果与某个“参考模型”预先给出的样本分组进行比较，即衡量聚类结果与预先分组情况的差异。内部指标只关注聚类的内部结构，即衡量簇内结构是否紧密，以及簇间距离是否拉开等。聚类分析和降维92

聚类分析和降维93

聚类分析和降维94

聚类分析和降维95

聚类分析和降维96

聚类分析和降维97

聚类分析和降维98

聚类分析和降维99

聚类分析和降维100K均值聚类算法（K-meansclusteringalgorithm）一种迭代求解的划分聚类方法（目前聚类中应用最广泛）其具体步骤如下：将数据随机分为K组，并随机选取K个对象作为每组的初始簇中心（或随机产生K个初始簇中心）对每一个样本，计算其与所有簇中心之间的距离，将其分配给距离它最近的簇中心如果没有点发生分配结果的改变或预先达到设定的收敛条件，则结束；否则进入下一步根据当前簇的样本重新计算并更新簇中心返回第2步聚类分析和降维101K均值聚类算法K均值算法的基本思想是让每个样本点离本簇中心的距离要小于与其他的距离优化目标是使每个样本点到本簇中心的距离平方和(sumofsquareddistance,SSD)尽量小。聚类分析和降维102

聚类分析和降维103层次聚类（HierarchicalClustering）一种很直观的算法它是假设类别之间存在层次结构，试图在不同层次实现对数据的划分，或将样本聚到层次化的类中层次聚类可以有聚合式或自下而上（用的比较多）对小的类别进行聚合的凝聚法，也有聚类自上向下把大的类别进行分割的分裂法聚类分析和降维104

聚类分析和降维105层次聚类聚类分析和降维106采用单链接方式的AGNES算法的具体流程：

聚类分析和降维107密度聚类密度聚类是基于密度（某样本点给定邻域内的其它样本点的数量）进行分簇思想：当邻域的密度达到指定阈值时，就将邻域内的样本点合并到本簇内，如果本簇内所有样本点的邻域密度都达不到指定阈值，则本簇划分完毕，进行下一个簇的划分。典型的密度聚类算法：DBSCAN（Density-BasedSpatialClusteringofApplicationswithNoise）及其派生算法密度聚类对非凸簇数据很有效聚类分析和降维108非凸簇数据示例图降维降维就是降低数据的维度，是缓解维数灾难的一个重要途径通过某种数学变换将高维数据空间转变到一个低维“子空间”，在这个子空间中样本的密度大幅提高，距离计算更为容易。可以实现在减少不相关或冗余特征的同时，提高模型精确度和运行效率聚类分析和降维109降维数据降维就是寻找一个映射函数f:x→y，将高维向量x映射成低维向量

y按照是否有使用样本标签，将降维算法分为有监督降维（LDA算法），和无监督降维（主成分分析法）按照降维算法使用的映射函数，可以将算法分为线性降维与非线性降维（奇异值分解），而核主成分分析方法则属于典型的非线性降维方法非线性降维（主要介绍）通常通过核映射和流形学习两种方式实现聚类分析和降维110

聚类分析和降维111降维非线形降维——流形学习（ManifoldLearning）传统欧式空间的度量难以用于真实世界的非线性数据，需要对数据的分布引入新的假设流形学习是从高维采样数据中恢复低维流形结构，即找到高维空间中的低维流形，并求出相应的嵌入映射，以实现维数约简或者数据可视化。流形学习基于这样一种假设：若低维流形嵌入到高维空间中，则数据样本在高维空间的分布虽然看上去十分复杂，但局部上仍具有欧氏空间的性质，可以在局部建立降维映射关系，然后再设法将局部映射关系推广到全局。聚类分析和降维112降维流形学习——等度量映射（Isomap）最经典的非线性映射降维方法之一主要目标：对于给定的高维流形，欲找到其对应的低维嵌入，使得高维流形上数据点间的近邻结构在低维嵌入中得以保持Isomap在计算高维流形上数据点间距离时，不是采用传统的欧式距离，而是采用微分几何中的测地线距离（或称为曲线距离），并且找到了一种用实际输入数据估计其测地线距离的算法。为了描述图2.9中两点A和B之间的距离，虚线为欧式距离，实线为测地线距离聚类分析和降维113图2.9测地极限距离与高维欧氏距离降维流形学习——等度量映射（Isomap）计算流形中的测地线距离：使用流形在局部与欧氏空间同胚的性质，对每个点找出其近邻点，然后建立一个近邻连接图（近邻点之间有连接，非近邻点之间无连接）将计算两点之间的测地线距离转变为近邻连接图上两点之间的最短路径问题（可采用Dijkstra或Floyd算法）得到距离矩阵后，就可以通过MDS算法获得样本在低维空间中的坐标聚类分析和降维114

聚类分析和降维115

聚类分析和降维116

聚类分析和降维117

聚类分析和降维118119大数据计算方法Chap3：线性方程组求解与矩阵分解基于LU分解的方程求解条件数与舍入误差分析对称正定矩阵的Cholesky分解QR分解与正交三角化线性最小二乘问题解法Outline120基于LU分解的方程求解方程求解的高斯消去法部分主元LU分解的推导排列阵、上/下三角阵方程部分主元LU分解的程序有关Matlab命令121

求解线性方程组

122

高斯消去法123用高斯消去法求解单个线性方程组有时需要求解仅右端向量不同的多个方程组线性方程组求解器124常微分方程组的例子：用向后欧拉法(backwardEuler)求解线性电路瞬态仿真问题，假设时间步长为常数

t举例125,1

1F写出A,B固定不变反复用高斯消去法求解不同右端项的方程组效率低能否节省计算量、复用中间结果?LU分解126x1b1x2b2使用LU分解!bxyy只做1次反复做实际中，还要

考虑选主元部分主元(列主元)求解线性方程组的一个例子部分(列)主元的高斯消去过程

127

(第1个方程乘0.3,-0.5,加到下面的方程上)”乘子”:当前列其他系数除以主元，即-0.3,0.5互为相反数

部分(列选)主元的高斯消去过程(第1个方程乘0.3,-0.5,下加)消元时,乘子为-0.04128最后一个方程

L=

部分(列选)主元LU分解A矩阵A与L,U的关系：UPA=LU129M2P2M1A=U由乘子填充;选主元行交换称为部分主元LU分解行交换阵P2消去矩阵M1,M2M2”乘子”:第一次为-0.3,0.5；第二次为-0.04

排列阵P:由单位阵经一系列的行交换而得.每行(列)上有且仅有一个1，其余为0例:向量p表示单位阵各行的新顺序，简洁地表示排列阵I(p,:)就是排列阵P如何计算PA？求解系数为排列阵的方程>>B=A(p,:)排列阵与上(下)三角阵PA

~

对A作行重排AP~对A作列重排>>p=[4132];>>

I=eye(4);P=I(p,:)(P为正交阵)1301:n的一个重排%时间短,省内存I(:,p)中p代表列的顺序因此就是PTpermutationmatrixx(p(1))=b(1)x(p(2))=b(2)...>>x(p)=b;

系数矩阵为上(下)三角阵的方程矩阵的LU分解得:上三角阵U，单位下三角阵L解Ux=b的两种”回代”(backsubs)排列阵与上(下)三角阵x=zeros(n,1);forj=n:-1:1x(j)=b(j)/U(j,j);

i=1:j-1;b(i)=b(i)-x(j)*U(i,j)endx=zeros(n,1);fork=n:-1:1

j=k+1:n;x(k)=(b(k)-U(k,j)*x(j))/U(k,k)end从b中依次减去U某列的倍数用U的行与解出的部分x作内积131b;x(这样,b不改变)L=LU分解的算法132算法3.1不选主元的LU分解

乘子原地存储：lupp程序“原地存储”方式最后两行算出L和Ufunction[L,U,p]=lupp(A)[n,n]=size(A);p=(1:n)';fork=1:n-1[r,m]=max(abs(A(k:n,k)));m=m+k-1;

if(A(m,k)~=0)

if(m~=k)A([km],:)=A([mk],:);p([km])=p([mk]);

end

i=k+1:n;A(i,k)=A(i,k)/A(k,k);j=k+1:n;A(i,j)=A(i,j)-A(i,k)*A(k,j);

endend最大值,位置索引交换两行算乘子更新未消去部分(kijversion)133L=tril(A,-1)+eye(n);U=triu(A);部分主元LU分解程序尽量用向量/矩阵运算保证算法不中断防止小主元导致的误差传播例子:在一个只有3位十进制精度的机器上求解如果不选主元,直接用高斯消去法为什么必须选主元?

134对上例,如果采用列主元高斯消去法为什么必须选主元?135算法的稳定性不同!非常接近准确解!其他选主元策略全主元：从未消去部分选最大矩阵元素，然后交换行、列行主元：从当前行选最大元素，然后交换列带阈值选主元：选某个大于最大元素倍(如=0.2)的元素Matlab中相关命令lu:

部分主元LU分解基于LU分解求解方程Ax=b:\(backwardslash):

X=A\B求解线性方程组AX=B基于部分主元LU分解的算法:mybslash.m是”\”的简化版(习题3.6)部分主元LU分解程序136稀疏阵稠密阵Lapack(BLAS)UMFPACK/CHOLMODPA=LU

LUx=PbSolveLy=PbSolve

Ux=yO(n3)O(n2)O(n2)A(准确度与运行时间的tradeoff)条件数与舍入误差分析复习矩阵范数敏感性分析与矩阵条件数考虑计算误差的分析算法稳定性137范数与矩阵条件数(复习)

(曼哈顿范数)(“最大”范数)(欧氏距离)

>>x=(1:4)/5>>norm1=norm(x,1)>>norm2=norm(x)>>norminf=norm(x,inf)138范数与矩阵条件数(复习)

01

139范数与矩阵条件数(复习)

条件数矩阵

条件数140F-范数:

(不是诱导范数)01

矩阵条件数141线性变换对

单元圆的扭曲

(排列阵)

(在p范数下)

算法稳定性与解的误差142

J.H.Wilkinson矩阵计算的舍入误差分析(绝对值)

(增长因子)

算法稳定性与解的误差

交换行,消元(采用截断舍入)解出:若不知道准确解,计算残差

143

注意：残差或相对残差不总反映结果的准确度某种相对残差相对残差小是算法稳定的必要条件

算法稳定性与解的误差144cond~106

矩阵条件数

算法稳定性再次说明:算法稳定性与问题敏感性共同影响解的误差

Cholesky分解对称正定矩阵的分解分解算法与应用145Cholesky分解定理1.如果实对称阵A的前n-1个顺序主子式

0

(LDLT分解)：其中L为单位下三角阵，D为对角阵2.若矩阵A同时正定，则存在非奇异下三角阵L，若限定L

的对角元>0,则此分解唯一In2x2case,forexample,whichimplies

146Cholesky分解算法可根据LU分解算法的kij版本进行修改，考虑对称性；原地存储:仅使用A的下三角部分，L的结果将其覆盖1.Fork=1ton1’.2.Fori=k+1ton3.4.Forj=k+1toi5.6.End7.End8.End注意：L不是单位

下三角=AFork=1tonFori=k+1tonEndEnd第2种算法147(不用选主元)每步计算都求得L中的某个元素!对称正定矩阵的Cholesky分解算法计算量分析(稠密阵)：仅需要n3/6次乘除法和差不多的加减法，是LU分解的一半计算量对于对称正定矩阵，可证明：所有n个平方根运算的操作数都为正数，算法可行Cholesky算法是稳定的(对舍入误差不敏感)，不需要选主元Matlab命令为chol,缺省得到上三角矩阵(R'*R=A),

可用于检测矩阵正定性解Ax=b:>>R=chol(A)>>L=chol(A,‘lower’)148(仅使用A的上三角部分元素)A=RTR

RTRx=bSolve

RTy=bSolveRx=y复杂度:O(n3)O(n2)O(n2)应用实例

000

149QR分解与正交三角化QR分解Gram-Schmidt正交化过程Householder变换Givens变换150“六大矩阵分解”A=QRA为m

n矩阵(先假设mn)Q为正交阵，R为上(右)三角阵正交阵的特点？元素值;正交性;转置;乘积;保2、F范数精简形式:Q为正交阵的前n列(列正交阵)若m<n定理:

对实矩阵𝑨,一定存在QR分解;若𝑨列满秩,且要求𝑹的对角元都>0,精简QR分解唯一矩阵的QR分解151AAAQQA=QRRR>>[Q,R]=qr(A);>>[Q,R]=qr(A,0);(用Cholesky分解唯一性证明)基于Gram-Schmidt正交化的算法O

计算过程保证了{qi}中任意两元素正交152fork

=1tonend基于Gram-Schmidt正交化的算法ClassicalG-S

procedurefork

=1tonforj

=1tok-1endendModifiedG-S

procedure计算量:

~mn2用剩下的向量做投影，而不是用原始的,提高准确度得到精简QR分解若A不列满秩，算法中断正交性丧失问题fork

=1tonforj

=1tok-1endend相邻的qk与qk+1之间正交性更好重正交化：153

基于Householder变换的算法

S

154

基于Householder变换的算法

可能同符号数相减!

修改后无抵消155

(2m次乘法,4mflops)

验证：基于Householder变换的算法

156完整QR分解算法157

(稳定的算法！)flops

完整QR分解算法158

基于Givens变换的算法

将二阶Givens阵

嵌入n阶单位阵:

159

基于Givens变换的算法例:

用Givens旋转变换进行消元每个Givens旋转阵用参数c,s刻画每次旋转变换仅影响向量两个元素则:

解:先针对第1,3分量构造二阶旋转矩阵,

接着对第1,4分量旋转

,仅影响矩阵的两行

160基于Givens变换的算法做若干次Givens变换实现对稠密矩阵用Givens旋转来正交三角化,比基于Householder变换的方法计算量约多50%与Householder变换需存u向量不同，Givens旋转需存储c,s参数值，存储用量约是前者的两倍适合处理稀疏矩阵A(填入现象，通过列交换减少填入)A实现第1列消元实现第2列消元

(不影响第1列)...得到R实现AE=QR求解稀疏线性最小二乘回归161线性最小二乘问题解法概述法方程方法利用QR分解解线性最小二乘列不满秩的情况162

求解线性最小二乘Amn

奇异163Ax

b

被拟合函数值基函数采样点

求解线性最小二乘正交变换保2范数设

n行

取等号条件:164n行

Ax

b

不满秩情况

=Ax

b165QR分解近似0

矩阵R1近奇异怎么求某一个解？解不唯一!(Householder变换)(这种满秩可能由数值误差造成,将其看成秩为1更合理)Rankrevealing问题做正交变换时适当地交换列在未消去的子阵中找2-范数最大的列,交换到当前列若第k步消去后,第k+1列的2-范数为0,则后续列均为0,且

k=rank(A)(rankrevealing!)将QRCP用于求𝑨𝒙

b的基本解等价于求解列主元QR分解>>x=A\b;166设

设,可取

(QRCP)R对角元绝对值递减Matlab有关命令qr命令能处理稠密或稀疏矩阵,允许m<n[Q,R]=qr(A);完整的QR分解,Householder,Givens[Q,R]=qr(A,0);简化形式的QR分解，若m<n同上[Q,R,P]=qr(A);A稠密,P为排列阵,AP=QR,QRCP算法[Q,R,P]=qr(A,'vector');R=qr(A);

A稠密时，看R的上三角部分[Q,R,P]=qr(A);A稀疏时,排列阵P是为了减少R中fill-in[C,R]=qr(S,B);[C,R,P]=qr(S,B);解稀疏最小二乘Sx

B\命令求解线性最小二乘问题，列重排的Householder变换

列不满秩会报warning，得到一个基本解

polyfit,多项式回归：p=polyfit(x,y,n)167内部使用LAPACK包稠密阵P排列向量,A(:,P)=QR不生成Q,效率高,x=R\C[G,y]=

planerot(x);Q=orth(A)168大数据计算方法Chap4：大型稀疏线性方程组的求解稀疏矩阵基本概念稀疏矩阵的分解算法高效率的稀疏矩阵直接解法迭代解法与共轭梯度法广义最小残量法(GMRES)Outline169稀疏矩阵基本概念基本概念存储格式简单的处理算法170稀疏矩阵(sparsematrix)存在大量零元素的矩阵Wilkinson’sdefinition:“matricesthat

allowspecialtechniquestotakeadvantage

ofthelargenumberofzeroelements.”很多实际仿真问题(电路网络,偏微分方程...)矩阵都是稀疏的非零元数、稀疏度、稠密度带状矩阵(bandmatrix):基本概念>>density=nnz(A)/prod(size(A))>>sparsity=1-density

带宽:非零元到主对角线的最大距离>>[i,j]=find(A)>>bandwidth=max(abs(i-j))171(也称之为”半带宽”)非零元数目为O(n)(不处理零元素)spy看非零元分布情况存储格式一般稀疏矩阵的存储结构三元组(COO格式)非零元的排列顺序可任意;压缩稀疏行(CSR格式)按第1行,第2行,…顺序存非零元;压缩稀疏列

aa1253467891011row11222333445coa1253467891011corow13691112

方便做逐行逐行

执行的操作172冗余:可能连续出现相同行编号prow为各行第一个元素位置aa1364927510811row12324132435pcol14681012

(CSC格式)用二维数组存储的矩阵不是稀疏矩阵!Matlab中的稀疏矩阵采用压缩稀疏列(CSC)方式存储与稠密矩阵(二维数组)转换若阶数很大,稠密格式A无法存

可直接用sparse命令生成

spdiags:“带状”稀疏矩阵sprand,sprandn生成随机的稀疏阵>>whosNameSizeBytesClassAttributesA20x203200doubleS20x20528doublesparse基本的Matlab函数>>S=sparse(i,j,x,m,n)S=spdiags([abc],[-1,0,1],n,n)173>>S=sparse(A)>>A=full(S)(相当于用COO格式输入数据)生成的S满足whos看变量类型(相当于用DIA格式输入数据)带状矩阵的分解与矩阵向量乘基本原则通过不操作0减少计算量(与0的加减乘除都没有增加新信息)例如：A+B的浮点运算次数不超过LU/Chol等复杂的矩阵分解,应设计考虑稀疏矩阵的特殊算法三对角矩阵方程Ax=b

tridisolve.m,用3个一维数组存A1741.Fork=1ton-12.Fori=k+1ton3.4.Forj=k+1ton5.6.End7.End8.End主元乘子ak,k+2=0复杂度O(n)!min(nnz(A),nnz(B))LU分解带状矩阵的分解与矩阵向量乘带状矩阵的LU分解仍然保持原始带宽

很多应用问题中的矩阵A对角占优/对称正定，不需要选主元倘若需要选主元，带宽增大，最多使带宽达2

+11751.Fork=1ton-12.Fori=k+1tok+

-13.4.Forj=k+1tok+

-15.6.End7.End8.End~线性复杂度(

<<n时)

00乘法次数:带状矩阵的分解与矩阵向量乘矩阵与向量乘法的算法矩阵存储采用CSC或其变种C语言算法实现(y=Ax)176jrow[k]

x

yj

Avoid

matvecC(cscptr

mata,double*x,double*y){

intn=mata->n,j,k,*ki;

double*kv;

for(j=0;j<n;j++)y[j]=0;

for(j=0;j<n;j++){

jval=mata->value[j];

jrow=mata->rowi[j];for(k=0;k<mata->nzcount[j];k++)y[jrow[k]]+=jval[k]*x[j];}

return;}structcsc{intn;int*nzcount;int**rowi;double**value;}*cscptr//第j列非零元行号数组//第j列的非零元素数组复杂度O(nnz(A))!j稀疏矩阵的直接解法直接解法概述Lsolve算法稀疏矩阵的LU分解稀疏矩阵的Cholesky分解消去树与并行算法177178高斯消去过程中的“填入”1.满足何条件时(i,j)处会发生fill-in?2.若非零元分布对称,发生fill-in的位置也对称?XXXXXXXXXX填入给LU/Chol分解带来麻烦？!nz=4991nz=348902基本问题——填入(fill-in)稀疏方程的直接解法概述高斯消元法(LU/Chol分解)针对稀疏矩阵的改进如何利用稀疏性减少内存用量和计算时间？减少内存用量稀疏矩阵存储格式(CSR,CSC,…)考虑矩阵非零元分布的变化(fill-in)，精确分配存储空间减少计算量忽略与0元素有关的运算填入元(fill-in)使矩阵“不稀疏”(也影响内存量)矩阵重排序(reordering)、预排序(preordering)难点--如何同时考虑减少填入与数值选主元?179Fill-reductionnumericalpivotingLsolve算法

180

xjx,b

l:,j=(须反映计算xi的先后依赖关系)算x3依赖x1算x4依赖x2算x8依赖x3,x5......

Lsolve算法181

(1)(2)

记为X=ReachG(B)顶点j顶点iLsolve算法182有向边说明了xi计算的依赖关系Lsolve算法183遍历所有路径，标记经过的点有向边说明了xi计算的依赖关系Lsolve算法184遍历所有路径，标记经过的点遍历某顶点的出边，相当于遍历矩阵的某列非零元(矩阵采用CSC存储)有向边说明了xi计算的依赖关系Lsolve算法185注意:采用DFS进行后序遍历(栈),得到X的顺序正确(DAG的拓扑顺序)有向边说明了xi计算的依赖关系Lsolve算法1861234B={1}X={4}X={2,4}X={3,2,4}X={1,3,2,4}另一个例子DFS后序遍历得出的顺序是正确的(不违反xi之间的依赖关系)Lsolve算法187时间复杂度为O(flop);

flop可能比n和nnz(L)小得多j的邻点集合

(CSC存储结构)DFS后序遍历和栈稀疏LU分解算法左视LU分解算法(jki版本)188kij版本LU分解Right-looking

Left-looking被乘数乘子

189稀疏LU分解算法

只更新排列向量稀疏Cholesky分解算法190逐列计算的Cholesky分解算法

LLTAj由Left-lookingLU分解改造得来，只用、也只算下三角部分

利用Lsolve算法

或者，非零元位置预测可从循环中提出来

符号分解符号分解中的消去树191

顶点i顶点j顶点k

(j是i的双亲节点，节点n是树根)192符号分解中的消去树消去树的例子193符号分解中的消去树消去树的例子

etree_demo.m注意:对选主元的GPLU算法，由于交换行,不能用消去树194

195并行LU分解算法

快速计算符号分解反映了依赖关系的上限

(不受选主元影响)计算顺序反映依赖关系(比Lsolve算法中xi间依赖关系复杂)(串行)高效稳定的稀疏矩阵直接解法减少填入的矩阵重排稀疏对称正定阵的解法高效稳定的稀疏矩阵LU分解196集成电路供电线网仿真DC分析：求解电阻网络

电路中的节点电压一个IBM公司芯片设计的例子n=474,524nnz=2,020,882在Matlab中约占36MB内存“\”求解(SPD阵,内部用chol)

t

2s“R=chol(A)”:时间长,或内存

不够用!

?集成电路供电线网的一部分197?减少填入的重排—实例如何进行重排序？两个行/列对称重排的例子198Fill-ins交换1,2

行/列Nofill-in!(1,2,…,n-1,n)(n,2,…,n-1,1