Grassmann.jl A\b 3x faster than Julia's StaticArrays.jl

栏目: IT技术 · 发布时间: 4年前

内容简介:In this algebra, it’s possible to compute on a mesh of arbitrary 5 dimensionalAdditionally, inProgramming the

In this algebra, it’s possible to compute on a mesh of arbitrary 5 dimensional conformal geometric algebra simplices, which can be represented by a bundle of nested dyadic tensors.

julia> using Grassmann, StaticArrays; basis"+-+++"
(⟨+-+++⟩, v, v₁, v₂, v₃, v₄, v₅, v₁₂, v₁₃, v₁₄, v₁₅, v₂₃, v₂₄, v₂₅, v₃₄, v₃₅, v₄₅, v₁₂₃, v₁₂₄, v₁₂₅, v₁₃₄, v₁₃₅, v₁₄₅, v₂₃₄, v₂₃₅, v₂₄₅, v₃₄₅, v₁₂₃₄, v₁₂₃₅, v₁₂₄₅, v₁₃₄₅, v₂₃₄₅, v₁₂₃₄₅)

julia> value(rand(Chain{V,1,Chain{V,1}}))
5-element StaticArrays.SArray{Tuple{5},Chain{⟨+-+++⟩,1,  ,253} where 253 where   ,1,5} with indices SOneTo(5):
   -0.33459594357756073v₁ - 0.3920064467082769v₂ - 0.575474920388841v₃ + 0.6150287650825268v₄ - 0.7568209093000915v₅
  -0.7402635950699139v₁ - 0.9303076362461833v₂ + 0.9729806462365271v₃ - 0.8514563480551867v₄ + 0.09906887873006287v₅
  -0.7456570397821101v₁ - 0.6497560949330929v₂ + 0.4756585550844967v₃ - 0.31169948016530347v₄ - 0.9355928499340793v₅
   -0.4555014543082292v₁ + 0.712268225360094v₂ - 0.7500443783398549v₃ - 0.36349628003234713v₄ + 0.5005769037801056v₅
 -0.07402971220645727v₁ + 0.19911765119918146v₂ - 0.4980618129231722v₃ - 0.7728564279829087v₄ + 0.9165735719353756v₅

julia> A = Chain{V,1}(rand(SMatrix{5,5}))
(0.9244801277294266v₁ + 0.029444337884018346v₂ + 0.745495522394158v₃ + 0.6695874677964055v₄ + 0.4998003712198389v₅)v₁ + (0.5423877012973404v₁ + 0.30112324458605655v₂ + 0.9530587650033631v₃ + 0.2706004745612134v₄ + 0.37762612797501616v₅)v₂ + (0.7730171467954035v₁ + 0.019660709510785912v₂ + 0.39119534821037494v₃ + 0.9403026278575068v₄ + 0.07545094732793833v₅)v₃ + (0.7184128110093908v₁ + 0.6295740775044767v₂ + 0.5179035493253021v₃ + 0.039081667453648716v₄ + 0.3719284661613145v₅)v₄ + (0.5033657705978616v₁ + 0.41183905359914386v₂ + 0.7761548051732969v₃ + 0.07635587137916744v₄ + 0.5582934197259402v₅)v₅

Additionally, in Grassmann.jl we prefer the nested usage of pure ChainBundle parametric types for large re-usable global cell geometries, from which local dyadics can be selected.

Programming the A\b method is straight forward with some Julia language metaprogramming and Grassmann.jl by first instantiating some Cramer symbols

Base.@pure function Grassmann.Cramer(N::Int)
    x,y = SVector{N}([Symbol(:x,i) for i ∈ 1:N]),SVector{N}([Symbol(:y,i) for i ∈ 1:N])
    xy = [:(($(x[1+i]),$(y[1+i])) = ($(x[i])∧t[$(1+i)],t[end-$i]∧$(y[i]))) for i ∈ 1:N-1]
    return x,y,xy
end

These are exterior product variants of the Cramer determinant symbols ( N! times N -simplex hypervolumes), which can be combined to directly solve a linear system:

@generated function Base.:\(t::Chain{V,1,<:Chain{V,1}},v::Chain{V,1}) where V
    N = ndims(V)-1 # paste this into the REPL for faster eval
    x,y,xy = Grassmann.Cramer(N)
    mid = [:($(x[i])∧v∧$(y[end-i])) for i ∈ 1:N-1]
    out = Expr(:call,:SVector,:(v∧$(y[end])),mid...,:($(x[end])∧v))
    return Expr(:block,:((x1,y1)=(t[1],t[end])),xy...,
        :(Chain{V,1}(getindex.($(Expr(:call,:./,out,:(t[1]∧$(y[end])))),1))))
end

Which results in the following highly efficient @generated code for solving the linear system,

(x1, y1) = (t[1], t[end])
(x2, y2) = (x1 ∧ t[2], t[end - 1] ∧ y1)
(x3, y3) = (x2 ∧ t[3], t[end - 2] ∧ y2)
(x4, y4) = (x3 ∧ t[4], t[end - 3] ∧ y3)
Chain{V, 1}(getindex.(SVector(v ∧ y4, (x1 ∧ v) ∧ y3, (x2 ∧ v) ∧ y2, (x3 ∧ v) ∧ y1, x4 ∧ v) ./ (t[1] ∧ y4), 1))

Benchmarks with that algebra indicate a 3x faster performance than SMatrix for applying A\b to bundles of dyadic elements.

julia> @btime $(rand(SMatrix{5,5},10000)).\Ref($(SVector(1,2,3,4,5)));
  2.588 ms (29496 allocations: 1.44 MiB)

julia> @btime $(Chain{V,1}.(rand(SMatrix{5,5},10000))).\$(v1+2v2+3v3+4v4+5v5);
  808.631 μs (2 allocations: 390.70 KiB)

julia> @btime $(SMatrix(A))\$(SVector(1,2,3,4,5))
  150.663 ns (0 allocations: 0 bytes)
5-element SArray{Tuple{5},Float64,1,5} with indices SOneTo(5):
 -4.783720495603508
  6.034887114999602
  1.017847212237964
  6.379374861538397
 -4.158116538111051

julia> @btime $A\$(v1+2v2+3v3+4v4+5v5)
  72.405 ns (0 allocations: 0 bytes)
-4.783720495603519v₁ + 6.034887114999605v₂ + 1.017847212237964v₃ + 6.379374861538393v₄ - 4.1581165381110505v₅

Such a solution is not only more efficient than Julia’s StaticArrays.jl method for SMatrix , but is also useful to minimize allocations in Grassmann.jl finite element assembly.

Similarly, the Cramer symbols can also be manipulated to invert the linear system or for determining whether a point is within a simplex.

julia> using Grassmann; basis"3"
(⟨+++⟩, v, v₁, v₂, v₃, v₁₂, v₁₃, v₂₃, v₁₂₃)

julia> T = Chain{V,1}(Chain(v1),v1+v2,v1+v3)
(1v₁ + 0v₂ + 0v₃)v₁ + (1v₁ + 1v₂ + 0v₃)v₂ + (1v₁ + 0v₂ + 1v₃)v₃

julia> barycenter(T) ∈ T, (v1+v2+v3) ∈ T
(true, false)

There are multiple equivalent ways of computing the same results using the and : dyadic products.

julia> T\barycenter(T) == inv(T)⋅barycenter(T)
true

julia> sqrt(T:T) == norm(SMatrix(T))
true

The following Makie.jl streamplot was generated with the Grassmann.Cramer method and interpolated from Nedelec edges of a Maxwell finite element solution.

More info about these examples is at https://grassmann.crucialflow.com/dev/tutorials/dyadic-tensors

Hermann Grassmann was the inventor of linear algebra as we know it today.


以上所述就是小编给大家介绍的《Grassmann.jl A\b 3x faster than Julia's StaticArrays.jl》,希望对大家有所帮助,如果大家有任何疑问请给我留言,小编会及时回复大家的。在此也非常感谢大家对 码农网 的支持!

查看所有标签

猜你喜欢:

本站部分资源来源于网络,本站转载出于传递更多信息之目的,版权归原作者或者来源机构所有,如转载稿涉及版权问题,请联系我们

思考,快与慢

思考,快与慢

[美] 丹尼尔·卡尼曼 / 胡晓姣、李爱民、何梦莹 / 中信出版社 / 2012-7 / 69.00元

《纽约时报》2011年度十大好书 新书上市,连续20多周蝉联亚马逊、《纽约时报》畅销书排行榜前20名,上市至今超过7个月,横扫全球各大畅销书排行榜,稳居亚马逊总榜前50名 《经济学人》、《华尔街日报》、《卫报》、《纽约时报》、《金融时报》、《商业周刊》、《华盛顿邮报》、等国外权威媒体,《三联生活周刊》、《商学院》、《东方早报》等国内知名媒体争相报道,国内外读者好评如潮 人类究竟有......一起来看看 《思考,快与慢》 这本书的介绍吧!

JS 压缩/解压工具
JS 压缩/解压工具

在线压缩/解压 JS 代码

Base64 编码/解码
Base64 编码/解码

Base64 编码/解码

html转js在线工具
html转js在线工具

html转js在线工具