diff --git a/notebooks/circular.jl b/notebooks/circular.jl index 9a85919..8b2777b 100644 --- a/notebooks/circular.jl +++ b/notebooks/circular.jl @@ -15,11 +15,14 @@ macro bind(def, element) end # ╔═╡ 6486e004-c95b-4fcd-936e-6508a0e0283c -using LinearAlgebra, StatsBase, PlutoUI +using LinearAlgebra, StatsBase, PlutoUI, GenericLinearAlgebra # ╔═╡ fab8681e-0be3-4d97-a0cb-2ad57472282a using Plots +# ╔═╡ 5fcb9846-0a0d-4517-8eff-90f580172b82 +using Quaternions + # ╔═╡ 5455f492-7c41-11ef-3315-19809ebaffc3 function circular(n,β) @@ -41,6 +44,9 @@ UpperHessenberg(z) end +# ╔═╡ db8c16a0-79df-4cc0-9c69-a6bcfcd23a1c +circular(3,2) + # ╔═╡ f9f00c9d-633a-4734-838f-3ffdcef6efe0 @bind β Slider([.5,1,4,10,100,1000,1e10]; default=1, show_value=true) @@ -68,53 +74,163 @@ function circular2(n) end # ╔═╡ 148317c1-e48f-4d39-833c-b51d2c43e5fd -function circular1(n) +function circular1(n) #β=1 A = randn(n,n) .+ im * randn(n,n) U = qr(A).Q * transpose( Matrix(qr(A).Q )) end -# ╔═╡ e14122da-9861-401c-bb56-3693b25e6e15 -circular2(5) - # ╔═╡ 5cd13d1b-9bb4-4a19-9e6a-07e214572639 +# ╠═╡ disabled = true +#=╠═╡ let n = 5 - t = 50_000 + t = 50 m = [ smallabs(circular2(n)) for i=1:t] display((mean(m),var(m))) stephist(m,normalize=true,label="dense") mh = [ smallabs(circular(n,2)) for i=1:t] display((mean(mh),var(mh))) - stephist!(mh,normalize=true,label="Hessenberg") - + title!("β=2 n=$n t=$t") + stephist!(mh,normalize=true,label="Hessenberg") end + ╠═╡ =# # ╔═╡ 63bf73bf-b312-49cc-a909-c043bfc80fd1 circular1(5) # ╔═╡ 82ce98e9-37a9-4023-ac62-5c092a8fa932 +# ╠═╡ disabled = true +#=╠═╡ let - n = 5 - t = 50_000 + n = 10 + t = 50 m = [ smallabs(circular1(n)) for i=1:t] display((mean(m),var(m))) stephist(m,normalize=true,label="dense") mh = [ smallabs(circular(n,1)) for i=1:t] display((mean(mh),var(mh))) + title!("β=1 n=$n t=$t") + stephist!(mh,normalize=true,label="Hessenberg") +end + ╠═╡ =# + +# ╔═╡ 686fae20-6ed7-444c-9085-0f97d8432303 +function circularm(n,β) + e = eigvals(circular(n,β)) + m2 = sum(e.^2) + m11 = (sum(e)^2-m2)/2 + (m2,m11) +end + +# ╔═╡ 5ede4dd3-c432-4108-b22d-ab72200f91f1 + + +# ╔═╡ aaa6835c-7d8f-4208-ab5d-0b06de61dbb8 +let + n = 3 + t = 10000 + β = 5 #-4/5 + # β = 1 # -1 + # β = 2 # -1 + v = [ circularm(n,β) for i=1:t] + # m2 = mean(first.(v)) + # m11 = mean(last.(v)) + # mean( first.(v) .* conj.(last.(v))) + # 3/mean(abs.(last.(v)).^2) - 1 # 3/|m11|²-1 = β + + mean(first.(v).*conj(last.(v))) , -6β/( (β+1)*(β+2)) +end + +# ╔═╡ dd7e4afa-323e-4287-a8a0-f1448fb669ac +5/7 + +# ╔═╡ fb3ec55c-db9f-4897-b4d0-3a07a3ba859c +# ╠═╡ disabled = true +#=╠═╡ +let + n = 5 + β = 50000 + t = 50 + mh = [ smallabs(circular(n,β)) for i=1:t] * 180/π + display((mean(mh),var(mh))) + stephist(mh,normalize=true,label="Hessenberg") + title!("β=$β n=$n t=$t") +end + ╠═╡ =# + +# ╔═╡ 37236855-5f7d-4d4a-8046-116a51dbae7f +function circular4(n) #β=4 + A = randn(2n,2n) + im * randn(2n, 2n) + U = Matrix(qr(A).Q).* cis.(2π*rand(2n)) + Z = kron( [0 -1;1 0], Matrix(I,n,n) ) + -Z * transpose(U) * Z * U +end + +# ╔═╡ 2443d72e-42f4-4f66-8170-04633f9b9031 +M = circular4(2) + +# ╔═╡ 662e1079-13d9-4c26-b359-96556bbf75ee + ℍ2ℂ(q::Quaternion) = [q.s+im*q.v1 q.v2+im*q.v3 + -q.v2+im*q.v3 q.s-im*q.v1 ] + +# ╔═╡ b11d81bd-b6d7-4325-9012-c65efb859d02 +ℍ2ℂ(M::Matrix) = hvcat(size(M,1),ℍ2ℂ.(M)...) + +# ╔═╡ a173bcd2-c02a-4897-9fbe-b7bde52f354b +negatej(q::Quaternion) = Quaternion(q.s,q.v1,-q.v2,q.v3) + +# ╔═╡ 61e9b1b5-95c5-4974-83b6-53387e93ded6 +function jdual(M) + negatej.(transpose(M)) +end + +# ╔═╡ b58a06f2-ad63-4742-abf8-bf491842b1f2 +function circular4q(n) #β=4 + A = quat.( randn(n,n), randn(n,n), randn(n,n), randn(n,n)) + U = quat(0,0,1,0) * Matrix(qr(A).Q) * quat(0,0,1,0) * jdual(Matrix(qr(A).Q)) #<-- this might be right +end + +# ╔═╡ 5e4d9a97-db38-445e-8982-7813b9656546 +let + n = 5 + t = 5000 + m = [ smallabs(circular4(n)) for i=1:t] + + stephist(m,normalize=true,label="dense") + mh = [ smallabs(circular(n,4)) for i=1:t] + mq = [ smallabs(ℍ2ℂ(circular4q(n))) for i=1:t] + title!("β=4 n=$n t=$t") stephist!(mh,normalize=true,label="Hessenberg") + stephist!(mq,normalize=true,label="Quaternionic Way") end +# ╔═╡ d426a597-d198-4d16-9ca1-4600c84d968f +Z = circular4q(3) + +# ╔═╡ b17edbda-3a23-48a1-8413-2a5350a29e86 +eigvals(ℍ2ℂ(Z)) + +# ╔═╡ 228f2a0c-2a52-432d-9dda-45f22dda4a4a +A = [quat( [rand(1:5) for i=1:4]...) for i=1:2,j=1:2] + +# ╔═╡ bc6bb958-7b18-4324-bc1d-644d543f9842 +jdual(A) + # ╔═╡ 00000000-0000-0000-0000-000000000001 PLUTO_PROJECT_TOML_CONTENTS = """ [deps] +GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8" +Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" [compat] +GenericLinearAlgebra = "~0.3.13" Plots = "~1.40.8" PlutoUI = "~0.7.60" +Quaternions = "~0.7.5" StatsBase = "~0.34.3" """ @@ -124,7 +240,7 @@ PLUTO_MANIFEST_TOML_CONTENTS = """ julia_version = "1.10.4" manifest_format = "2.0" -project_hash = "7b8203dbbe986d66cb4ddfe6db1b6ccaff902edb" +project_hash = "d5477210c2247950b3cde5a6471f10197fb68f21" [[deps.AbstractPlutoDingetjes]] deps = ["Pkg"] @@ -339,6 +455,12 @@ git-tree-sha1 = "a8863b69c2a0859f2c2c87ebdc4c6712e88bdf0d" uuid = "d2c73de3-f751-5644-a686-071e5b155ba9" version = "0.73.7+0" +[[deps.GenericLinearAlgebra]] +deps = ["LinearAlgebra", "Printf", "Random", "libblastrampoline_jll"] +git-tree-sha1 = "f47136cac29a9b7a8a88dbce1195394978091edb" +uuid = "14197337-ba66-59df-a3e3-ca00e7dcff7a" +version = "0.3.13" + [[deps.Gettext_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"] git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046" @@ -778,6 +900,12 @@ git-tree-sha1 = "729927532d48cf79f49070341e1d918a65aba6b0" uuid = "e99dba38-086e-5de3-a5b1-6e4c66e897c3" version = "6.7.1+1" +[[deps.Quaternions]] +deps = ["LinearAlgebra", "Random", "RealDot"] +git-tree-sha1 = "9a46862d248ea548e340e30e2894118749dc7f51" +uuid = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0" +version = "0.7.5" + [[deps.REPL]] deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" @@ -786,6 +914,12 @@ uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" deps = ["SHA"] uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +[[deps.RealDot]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "9f0a1b71baaf7650f4fa8a1d168c7fb6ee41f0c9" +uuid = "c1ae055f-0cd5-4b69-90a6-9a35b1a98df9" +version = "0.1.0" + [[deps.RecipesBase]] deps = ["PrecompileTools"] git-tree-sha1 = "5c3d09cc4f31f5fc6af001c250bf1278733100ff" @@ -1249,15 +1383,33 @@ version = "1.4.1+1" # ╠═6486e004-c95b-4fcd-936e-6508a0e0283c # ╠═fab8681e-0be3-4d97-a0cb-2ad57472282a # ╠═5455f492-7c41-11ef-3315-19809ebaffc3 +# ╠═db8c16a0-79df-4cc0-9c69-a6bcfcd23a1c # ╠═f9f00c9d-633a-4734-838f-3ffdcef6efe0 # ╠═e4418523-234c-4faa-a29b-43484013b322 # ╠═63cd9672-d3d2-4b2a-8076-f98b60dc11e6 # ╠═02db4122-8e3c-477c-b1ef-d9c236787910 # ╠═14b4a088-ff74-43f3-b301-4057593c0447 # ╠═148317c1-e48f-4d39-833c-b51d2c43e5fd -# ╠═e14122da-9861-401c-bb56-3693b25e6e15 +# ╠═2443d72e-42f4-4f66-8170-04633f9b9031 # ╠═5cd13d1b-9bb4-4a19-9e6a-07e214572639 # ╠═63bf73bf-b312-49cc-a909-c043bfc80fd1 # ╠═82ce98e9-37a9-4023-ac62-5c092a8fa932 +# ╠═5e4d9a97-db38-445e-8982-7813b9656546 +# ╠═d426a597-d198-4d16-9ca1-4600c84d968f +# ╠═b17edbda-3a23-48a1-8413-2a5350a29e86 +# ╠═686fae20-6ed7-444c-9085-0f97d8432303 +# ╠═5ede4dd3-c432-4108-b22d-ab72200f91f1 +# ╠═aaa6835c-7d8f-4208-ab5d-0b06de61dbb8 +# ╠═dd7e4afa-323e-4287-a8a0-f1448fb669ac +# ╠═5fcb9846-0a0d-4517-8eff-90f580172b82 +# ╠═fb3ec55c-db9f-4897-b4d0-3a07a3ba859c +# ╠═37236855-5f7d-4d4a-8046-116a51dbae7f +# ╠═662e1079-13d9-4c26-b359-96556bbf75ee +# ╠═b11d81bd-b6d7-4325-9012-c65efb859d02 +# ╠═a173bcd2-c02a-4897-9fbe-b7bde52f354b +# ╠═61e9b1b5-95c5-4974-83b6-53387e93ded6 +# ╠═b58a06f2-ad63-4742-abf8-bf491842b1f2 +# ╠═228f2a0c-2a52-432d-9dda-45f22dda4a4a +# ╠═bc6bb958-7b18-4324-bc1d-644d543f9842 # ╟─00000000-0000-0000-0000-000000000001 # ╟─00000000-0000-0000-0000-000000000002