Julia数据科学一览：线性代数和图像处理(1)

前言

 

去年夏天的时候，写过一篇《为什么说Julia比Python高出一个境界》，我在里面着重介绍了Julia，作为一种专门为数据科学而生的语言，其在速度，效率，优雅，易学等方面优势明显而且前程远大。本系列，先不进入语言细节，而是采取专题形式，首先向读者展示一下Julia在数据科学一些重要方面的实践内容。如果你有Python或者R语言，以及一些数据分析和数据挖掘的基础，可以相对轻松阅读以下相关内容，并能够体会到Julia语言在这其中体现的价值。在数据科学一览之后，我会以此为基础，来带大家回溯一下这些专题中的具体细节和重点内容.如果必要的话，之后还会做一个Julia语言的入门教程。

 

本篇用的julia版本为1.53，是目前的最新稳定版，每次升级版本时，我会在篇前说明。

 

线性代数

 

Julia语言的主打之一，便是科学计算，而线性代数又是科学计算的重要基础之一，我们先大致看一下julia在线性代数方面的表现。当然线性代数内容庞大，这里仅仅略举几个栗子，后续有空的话，可以开一个专门讲解线性代数的课程。

在计算机中，矩阵的存储方式如下图所示

先把可能用到的包载入

using LinearAlgebra #线性代数using SparseArrays  #稀疏数组using Images  #图像处理using MAT #  读取MATLAB文件

 

基本矩阵和线性方程组

 

创建一个随机矩阵，求其转置，再讲二者相乘，得到一个对称阵，其中代码末尾的分号意味着结果会被隐去，所以第3行没有结果显示，所显示的第4行的结果

A = rand(10,10); # created a random matrix of size 10-by-10Atranspose = A' # matrix transposeA = A*Atranspose; # matrix multiplicationA

结果显示如下，是一个float类型的二维矩阵

 

这里用一个宏@show来显示一下最终的矩阵A结果中，对角元素的一个判断结果，此时分号不会隐藏结果

@show A[11] == A[1,2];

 

 

解一个常规的线性方程组，如下，生成系数向量b，第2行是给出求解公式，第3行使用norm函数进行实际求解，并显示

b = rand(10); #created a random vector of size 10x = A\b; #x is the solutions to the linear system Ax=b@show norm(A*x-b);

结果如下：

 

需要说一下，以上的A是矩阵类型(Matrix)，b是向量(Vector)，转置是所谓伴随类型(Adjoint)，从中可见Julia的强类型特色，用typeof再看一下各自的类型​​​​​​​

@show typeof(A)@show typeof(b)@show typeof(rand(1,10))@show typeof(Atranspose);

 

有关adjoint类型，可以用？查一下，有关查询文档，跟Python，R类似。

?adjoint

 

 

使用parent方法，可以复原其转置前的样子，显然这是跑的一次操作，而不是数据本身

Atranspose.parent

 

 

如果想把转置后的对象copy过来，用copy函数即可​​​​​​​

# To actually copy the matrix:B = copy(Atranspose)

结果如下，从这几个矩阵的大小来看，能够体会到一点Julia的设计理念吧

 

再说一点，从数值计算角度，在求线性方程组的时候，推荐用 \ 表示除法，而求矩阵的逆存在数值计算上不安全的因素，不推荐使用。

 

 

矩阵分解

 

以下几种分解矩阵的方式，在实用角度，都为了应付不同的线性方程组求解问题

LU分解，是把一个矩阵分解成一个下三角阵和上三角阵乘积的形式

 

调用lu函数，代码如下：

luA = lu(A)

 

结果如下：

 

QR分解：实际上对LU分解的进一步正交化，结果是得到正交阵和上三角阵

可以调用qr函数

qrA = qr(A)

 

结果如下：

Cholesky分解，该方法仅针对正定矩阵而言，结果为一个下三角阵与其转置的乘积

按照常规步骤，先来判断是否为正定阵

isposdef(A)

结果如下，满足条件，否者再找其他

 

这里调用相应函数

cholA = cholesky(A)

 

结果如下

 

 

类型计算

 

我们看一个专门做分解的函数，在参数中没有做任何分解类型的规定

factorize(A)

 

结果如下，是一个Cholesky分解

 

原因何在呢，我们看下这个函数的文档

?factorize

 

从下面的说明可以知道，该函数是根据不同的输入类型，采取不同的返回值求法，这是Julia语言计算高效的主要原因之一。

以下介绍下对角阵和单位阵的创建

 

对角阵创建​​​​​​​

# convert(Diagonal{Int64,Array{Int64,1}},diagm([1,2,3]))Diagonal([1,2,3])

 

单位阵创建

I(3)

 

 

可见都是很容易的，细节的地方，后续详谈

 

稀疏矩阵

 

由于稀疏矩阵本身的稀薄疏离的性质，其在计算中以压缩形式存储，这里载入相关的包，创建一个随机的稀疏矩阵​​​​​​​

using SparseArraysS = sprand(5,5,2/5)

结果如下

​​​​​​​

 

有关该函数的用法，以上代码的前两项是对应稀疏矩阵的行列，第三个参数是非零元素的比例。这里看下矩阵视角下的结果：

 

 

 

结语

 

本篇内容相对抽象和枯燥一点，但没有办法，这个是涉及线性代数不可避免的，下一章，我们会以本篇的知识和代码内容为基础，介绍一些图像处理的东西。或许有了些图，会读起来相对轻松一些。

 

 

主要参考资料

https://github.com/JuliaAcademy/DataScience

Data Science with Julia ，Paul D. McNicholas, Peter A. Tait，2019

Julia Quick Syntax Reference: A Pocket Guide for Data Science Programming

Antonello Lobianco，2019

Julia for Data Science，Zacharias Voulgaris，2016

Mastering Julia: Develop your analytical and programming skills further in Julia to solve complex data processing problems，Malcolm Sherrington，2015

​​​​​​​