julia系列4：数学计算

IE06

已于 2022-08-01 17:02:06 修改

阅读量1.6k

点赞数 1

分类专栏： Julia系列文章标签： julia python 机器学习

于 2020-08-12 08:12:28 首次发布

本文链接：https://blog.csdn.net/kittyzc/article/details/107140644

版权

Julia系列专栏收录该内容

17 篇文章 8 订阅

订阅专栏

1. 基础

1.1 Array

在julia中，vector和matrix都是array：

julia> Vector
Array{T,1} where T
 
julia> Matrix
Array{T,2} where T

我们还可以显式地声明下面地矩阵：Diagonal(对角矩阵)，Symmetric(对称矩阵)，Hermitian(埃尔米特矩阵)
Tridiagonal(三对角矩阵)，SymTridiagonal(对称三角矩阵)[SymTridiagonal只能作用于对称矩阵]

n = 6
A = randn(n,n)
B = Symmetric(A)
B = Hermitian(A)
B = Tridiagonal(A)
B = SymTridiagonal(A+A')

注意julia中的slice是copy操作，引用操作需要使用view。

z = view(v,3:5)
# or use macro
z = @views v[3:5]

1.2 Dataframe

# 新建一个DataFrame并增加4列内容

using DataFrames
df1 = DataFrame()
df1[:clo1] = Array([1.0,2.0,3.0])
df1[:clo2] = Array([4.0,5.0,6.0])
df1[:clo3] = Array([7.0,8.0,9.0])
df1[:ID] = Array(['a','b','c'])
show(df1)
>>3×4 DataFrame
│ Row │ clo1    │ clo2    │ clo3    │ ID   │
│     │ Float64 │ Float64 │ Float64 │ Char │
├─────┼─────────┼─────────┼─────────┼──────┤
│ 1   │ 1.0     │ 4.0     │ 7.0     │ 'a'  │
│ 2   │ 2.0     │ 5.0     │ 8.0     │ 'b'  │
│ 3   │ 3.0     │ 6.0     │ 9.0     │ 'c'  │

# 如果没有指定列名，则默认是x1,x2...

df2 = DataFrame(rand(5,6))
# 在DataFrame定义时直接指定内容

df3 = DataFrame([collect(1:3),collect(4:6)], [:A, :B])
>>    A   B
  Int64    Int64
1    1   4
2    2   5
3    3   6

在这里插入图片描述

RDatasets是Julia中的一个数据集，里面包含了很多可以学习和验证的数据，其中就包括iris数据集。

下面是直接读取csv、xlxs、txt文件的方法

在这里插入图片描述

下面是一些基础操作：

注意最后三者的差别。
下面是修改数据的方法：

1.3 时间序列

在这里插入图片描述

2. 性能和稳定性

2.1 简单说明

首先要关注一下特殊符号，例如eltype函数：

function eltype(::Type{<:AbstractDict{K,V}}) where {K,V}
    if @isdefined(K)
        if @isdefined(V)
            return Pair{K,V}
        else
            return Pair{K}
        end
    elseif @isdefined(V)
        return Pair{k,V} where k
    else
        return Pair
    end
end

其中：
1.{}用来表示类型的参数化（泛型），比如说Vector{Int}表示一个整数的向量，在Java（以及C++）中用<>表示，例如List
2.[]表示数组
3.()可以表示元组（和Haskell,Python一致），如（1，2，3），也可以表示函数调用，如sum([1,2,3])
3.A<:B表示类型子类化（A是B的子类，a::B表示a的类型是B）
4.where引导类型变量，例如Array{T,3} where T表示一个3维数组，数组元素的类型为某个待定的T，可以提高性能
5.一个点 . ，如A.a，表示A有个属性为a；. 还可以用来表示broadcast，可以理解为就是map，例如说，你有一个矩阵A，sin.(A)表示将A中元素分别求正弦后得到的新矩阵，等价于map(sin,A)
6.@表示宏调用（Julia本来也可以不用@标记宏，但是为了大家方便阅读，最后还是强制要求用@标记）
7.Julia中有生成函数(一种特殊函数），本质上是一个返回表达式的函数，然后表达式会插入函数所在处执行了。这个东西创造出来是为了利用LLVM的多层次编译优势（运行时动态生成优化代码），主要用于优化（而且被大量使用）例如在文档中的例子：

julia> @generated function bar(x)
           if x <: Integer
               return :(x ^ 2)
           else
               return :(x)
           end
       end
bar (generic function with 1 method)
 
julia> bar(4)
16
 
julia> bar("baz")
"baz"

在Julia中，值只有具体类型，不可能有抽象类型，只有抽象类型能被子类化，具体类型不可以，具体说来，就是上述类型树中，只有叶子节点才能被构造出来，其他的节点是不可构造的。这是很有趣的一个地方，这表明，Julia的tag系统全都是具体类型组成的，那抽象类型有什么用呢？抽象类型用来确定子类化关系，以做多重派发以实现ad hoc多态（并且只有这个用处）

避免全局变量，把全局变量声明为常量可以巨大的提升性能。如果必须要声明全局变量，可以在使用它的地方标注他们的类型来优化效率。
任何注重性能或者需要测试性能的代码都应该被放置在函数之中。
可以通过@code来查看用于code generation的宏
下面是测试例子：

global x = rand(10000)

function lp1() # 直接动用全局变量
    s = 0.0
    for i in x
        s += i
    end
end

function lp2() # 指定全局变量类型
    s = 0.0
    for i in x::Vector{Float64}
        s += i
    end
end

function lp3(x) # 将变量以参数形式传入
    s = 0.0
    for i in x
        s += i
    end
end

在这里插入图片描述
指定类型是速度最快的，原理如下：
当我们定义一个函数时，如果函数参数的类型是固定的，比如是一个Int64的Array[1,2,3,4],那他们在内存中会连续存放；
但如果函数参数的类型是Any，那么内存中连续存放只是他们的“指针”,会指向其实际的位置。这样一来，数据存取就慢下来了。计算concrete类型会比计算abstract类型要节省时间，我们可以使用@code_warntype来查看运行的函数中是否有abstract类型，对于有abstract类型的地方，会用红色的标出。
在这里插入图片描述
可以用如下函数来获得相同的数据类型：
zero(value)
eltype(array)
one(value)
similar(array)
向量化并不会提高Julia的运行速度，@simd 可以在运算支持被重新recorded时加速运算
输出预分配可以提高性能
避免不必要的Array，比如计算x,y,z的和时，使用x+y+z，不要用sum([x,y,z])
用div(x,y)代替trunc(x/y)，用fld(x,y)代替floor(x/y)，用cld(x,y)代替ceil(x/y)，有现成的函数就不要自己写计算过程

2.2 全局变量设置为const

注意下面两段代码的差别，性能差距将近100倍。const类型的值是可以修改的，只是不能修改类型。用@code_warntype 可以看出两者的差别， pow_array 返回值是 Any 类型，而 pow_array2 是 Float64 类型，可见 pow_array2 是 type-stable 的。

julia> p = 2;
julia> function pow_array(x::Vector{Float64})
         s = 0.0
         for y in x
           s = s + y^p
         end
         return s
       end
pow_array (generic function with 1 method)
 
julia> t = rand(100000);
 
julia> @benchmark pow_array(t)
BenchmarkTools.Trial:
  memory estimate:  4.58 MiB
  allocs estimate:  300000
  --------------
  minimum time:     7.385 ms (0.00% GC)
  median time:      8.052 ms (0.00% GC)
  mean time:        8.261 ms (2.76% GC)
  maximum time:     50.044 ms (85.05% GC)
  --------------
  samples:          604
  evals/sample:     1
julia> const p2 = 2
2

julia> function pow_array2(x::Vector{Float64})
         s = 0.0
         for y in x
           s = s + y^p2
         end
         return s
       end
pow_array2 (generic function with 1 method)
 
julia> @benchmark pow_array2(t)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     81.324 μs (0.00% GC)
  median time:      83.629 μs (0.00% GC)
  mean time:        87.973 μs (0.00% GC)
  maximum time:     185.029 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1

2.3 fastMath

注意下面两个函数的差别。@fastmath宏会放宽一些限制，包括检查 NaN 或 Infinity 等，类似于 GCC 中的-ffast-math编译选项。主要是把普通的加减乘除换成了Base.FastMath中的方法。差异还是蛮大的，差不多 8 倍差距！

julia> function sum_diff(x)
             n = length(x); d = 1/(n-1)
             s = zero(eltype(x))
             s = s +  (x[2] - x[1]) / d
             for i = 2:length(x)-1
                 s =  s + (x[i+1] - x[i+1]) / (2*d)
             end
             s = s + (x[n] - x[n-1])/d
         end
sum_diff (generic function with 1 method)
 
julia> function sum_diff_fast(x)
              n=length(x); d = 1/(n-1)
              s = zero(eltype(x))
              @fastmath s = s +  (x[2] - x[1]) / d
              @fastmath for i = 2:n-1
                  s =  s + (x[i+1] - x[i+1]) / (2*d)
              end
              @fastmath s = s + (x[n] - x[n-1])/d
          end
sum_diff_fast (generic function with 1 method)
 
julia> t=rand(2000);
 
julia> sum_diff(t)
1124.071808538789
 
julia> sum_diff_fast(t)
1124.0718085387887
 
julia> using BenchmarkTools
 
julia> @benchmark sum_diff(t)
BenchmarkTools.Trial:
  memory estimate:  16 bytes
  allocs estimate:  1
  --------------
  minimum time:     4.447 μs (0.00% GC)
  median time:      4.504 μs (0.00% GC)
  mean time:        4.823 μs (0.00% GC)
  maximum time:     17.318 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     7
 
julia> @benchmark sum_diff_fast(t)
BenchmarkTools.Trial:
  memory estimate:  16 bytes
  allocs estimate:  1
  --------------
  minimum time:     574.951 ns (0.00% GC)
  median time:      580.831 ns (0.00% GC)
  mean time:        621.044 ns (1.04% GC)
  maximum time:     65.017 μs (99.06% GC)
  --------------
  samples:          10000
  evals/sample:     183

2.4 数值稳定性

一般是看输入的微小变化，会对输出产生什么影响。对于一个映射 $y = f (x)$ 的计算，考虑其输出结果对输入的敏感度，定义条件数condition为
在这里插入图片描述
如果映射可微，则

以 $A x = b$ 为例条件数为：

julia可以用linearAlgebra的cond函数，python可以使用numpy.linalq中的cond函数

3. 线性代数库

在开始讲解之前，我们先建立一个简单的线性系统：
A = rand(3,3)
x = fill(1, (3,1))
b = A * x

3.1. LU分解

首先看LU分解：PA = LU。此处P为置换矩阵，L是下三角矩阵（对角线元素为1），U是上三角矩阵。
LU分解的意义在于消除Ax=b中的 $A^{-1}$ 计算。Ly=b和Ux=y可以通过代入法进行求解。
LU分解的步骤如图：
在这里插入图片描述

LU分解中的P是为了解决数据问题。在处理数值计算问题时候，一定要尽可能避免大数和小数相加减的问题。因此，在高斯消元的时候，做一些行置换（pivoting），是的一列中最大的数被调换到对角线上进行消元。

# julia使用函数lufact可以获得矩阵A,L,U
Alu = lu(A) 

# 查看Alu的type：
typeof(Alu)
P = Alu.P
L = Alu.L
U = Alu.U

当A是对称正定矩阵时，LU分解变为Cholesky分解，即 $U=L^T$ ，并且一定是数值稳定的。

3.2 QR分解

QR分解主要用于基的Gram–Schmidt正交化。A = QR. 此处Q是酉矩阵/正交矩阵，R是上三角阵

Aqr =  qr(A)
Q = Aqr.Q
R = Aqr.R

3.3 特征值分解

所谓的特征方向，就是在这个方向上，矩阵变换等价于线性变换。
特征值分解可以用来做降维。特征值比较小的方向可以直接去除。SVD也是对特征值分解，但是和特征值分解不同，SVD分解并不要求被分解的矩阵为方阵，定义矩阵A的SVD为如下： $A=U\Sigma V^T$

Asym = copy(A)
AsymEig = eig(Asym)
# 特征值与特征向量可以通过下面的命令获得
AsymEig.values
AsymEig.vectors

这里要注意一下：inv(AsymEig)*Asym = I

3.4 求解方程

可以使用HomotopyContinuation，地址为https://www.juliahomotopycontinuation.org/
下面是个例子：
在这里插入图片描述

# load the package
using HomotopyContinuation
# declare variables x and y
@var x y
# define the polynomials
f₁ = (x^4 + y^4 - 1) * (x^2 + y^2 - 2) + x^5 * y
f₂ = x^2+2x*y^2 - 2y^2 - 1/2
F = System([f₁, f₂])
result = solve(F)

4. 概率相关

4.1 基本用法

在这里插入图片描述

4.2 综合例子

rand：uniform分布
StatsBase.countmap：统计列表中每个元素的个数，生成一个字典
bar可以直接绘制字典
下面是个综合例子：
在这里插入图片描述

IE06

关注

1
点赞
踩
5

收藏

觉得还不错? 一键收藏
0
评论
julia系列4：数学计算

我们还可以显式地声明下面地矩阵Diagonal(对角矩阵)，Symmetric(对称矩阵)，Hermitian(埃尔米特矩阵)Tridiagonal(三对角矩阵)，SymTridiagonal(对称三角矩阵)[SymTridiagonal只能作用于对称矩阵]注意julia中的slice是copy操作，引用操作需要使用view。.........
复制链接

扫一扫

专栏目录