{"id":10133,"date":"2018-02-15T22:21:44","date_gmt":"2018-02-15T21:21:44","guid":{"rendered":"http:\/\/djalil.chafai.net\/blog\/?p=10133"},"modified":"2018-02-22T11:00:29","modified_gmt":"2018-02-22T10:00:29","slug":"playing-a-bit-with-julia","status":"publish","type":"post","link":"https:\/\/djalil.chafai.net\/blog\/2018\/02\/15\/playing-a-bit-with-julia\/","title":{"rendered":"Playing a bit with Julia"},"content":{"rendered":"

\"The
\nThis short post is devoted to a couple of Julia<\/a> programs.<\/p>\n

Gaussian Unitary Ensemble.<\/strong>\u00a0It is the Boltzmann-Gibbs measure with density proportional to
\n$$
\n(x_1,\\ldots,x_N)\\in\\mathbb{R}^N\\mapsto\\mathrm{e}^{-N\\sum_{i=1}^N x_i^2}\\prod_{i<j}(x_i-x_j)^2.
\n$$
\nIt is the distribution of the eigenvalues of a random $N\\times N$ Hermitian matrix distributed according to the Gaussian probability measure with density proportional to
\n$$
\nH\\mapsto\\mathrm{e}^{-N\\mathrm{Tr}(H^2)}.
\n$$
\nIt is a famous exactly solvable model of mathematical physics. There is a nice formula for the mean empirical spectral distribution $\\mathbb{E}\\mu_N$ where $\\mu_N=\\frac{1}{N}\\sum_{i=1}^N\\delta_{x_i}$, namely
\n$$
\nx\\in\\mathrm{R}\\mapsto\\frac{\\mathrm{e}^{-\\frac{N}{2}x^2}}{\\sqrt{2\\pi N}}\\sum_{\\ell=0}^{N-1}P_\\ell^2(\\sqrt{N}x)
\n$$
\nwhere ${(P_\\ell)}_{\\ell\\geq0}$ are the Hermite polynomials which are the orthonormal polynomials for the standard Gaussian distribution $\\mathcal{N}(0,1)$. The computation of the Laplace transform and a subtle reasonning, see this paper<\/a>, reveal that it converges as $N\\to\\infty$ towards the the Wigner semicircle distribution with density with respect to the Lebesgue measure
\n$$
\nx\\in\\mathbb{R}\\mapsto\\frac{\\sqrt{4-x^2}}{2\\pi}\\mathbf{1}_{x\\in[-2,2]}.
\n$$
\nHere is a nice plot followed by the Julia code used to produce it.
\n\"\"<\/p>\n

function normalized_hermite_polynomials_numeric(n,x)\r\n    if size(x,2) != 1\r\n        error(\"Second argument must be column vector.\")\r\n    elseif n == 0\r\n        return ones(length(x),1)\r\n    elseif n == 1\r\n        return x\r\n    else\r\n        P = [ ones(length(x),1) x ]\r\n        for k = 1:n-1\r\n            P = [ P  x.*P[:,k+1]\/sqrt(k+1)-P[:,k]*sqrt(1\/(1+1\/k)) ]\r\n        end\r\n        return P # matrix with columns P_0(x),...,P_n(x)\r\n    end\r\nend #function\r\n\r\nfunction gue_numeric(n,x)\r\n    if size(x,2) != 1\r\n        error(\"Second argument must be column vector.\")\r\n    else\r\n        y = sqrt(n)*x\r\n        p = normalized_hermite_polynomials_numeric(n-1,y).^2\r\n        return exp(-y.^2\/2) .* sum(p,2) \/ sqrt(2*pi*n)\r\n    end\r\nend #function\r\n\r\n# Pkg.add(\"Plots\") # http:\/\/docs.juliaplots.org\/latest\/tutorial\/\r\nusing Plots\r\npyplot()\r\nN = 8\r\nX = collect(-3:.01:3)\r\nY = gue_numeric(N,X)\r\nfunction semicircle(x)\r\n    if (abs(x)>2) return 0 else return sqrt(4-x^2)\/(2*pi) end\r\nend # function\r\nY = [Y [semicircle(x) for x in X]]\r\nplot(X,Y,\r\n     title = @sprintf(\"GUE mean Empirical Spectral Distribution beta=2 N=%i\",N),\r\n     titlefont = font(\"Times\", 10),\r\n     label = [\"GUE mean ESD\" \"SemiCircle\"],\r\n     lw = 2,\r\n     linestyle = [:dash :solid],     \r\n     linecolor = [:darkblue :darkred],\r\n     aspect_ratio = 2*pi,\r\n     grid = false,\r\n     border = false)\r\ngui()\r\n#savefig(\"gue.svg\")\r\nsavefig(\"gue.png\")\r\n<\/pre>\n
Here is, just for fun, a way to produce symbolically Hermite polynomials.<\/p>\n
# Pkg.add(\"Polynomials\") # http:\/\/juliamath.github.io\/Polynomials.jl\/latest\/\r\nusing Polynomials\r\n\r\nfunction unnormalized_hermite_polynomial_symbolic(n) \r\n    if n == 0\r\n        return Poly([1])\r\n    elseif n == 1\r\n        return Poly([0,1])\r\n    else\r\n        P = [Poly([1]) Poly([0,1])]\r\n        for k = 1:n-1\r\n            P = [ P Poly([0,1])*P[k+1]-k*P[k] ]\r\n        end\r\n        return P # better than just P[n+1]\r\n    end\r\nend #function\r\n<\/pre>\n
Complex Ginibre Ensemble.<\/strong> It is the Boltzmann-Gibbs measure with density proportional to
\n$$
\n(x_1,\\ldots,x_N)\\in\\mathbb{C}^N\\mapsto\\mathrm{e}^{-N\\sum_{i=1}^N|x_i|^2}\\prod_{i<j}|x_i-x_j|^2.
\n$$
\nIt is the distribution of the eigenvalues of a random $N\\times N$ complex matrix distributed according to the Gaussian probability measure with density proportional to
\n$$
\nM\\mapsto\\mathrm{e}^{-N\\mathrm{Tr}(MM^*)}
\n$$
\nwhere $M^*=\\overline{M}^\\top$. Yet another exactly solvable model of mathematical physics, with a nice formula for the mean empirical spectral distribution $\\mathbb{E}\\mu_N$ where $\\mu_N=\\frac{1}{N}\\sum_{i=1}^N\\delta_{x_i}$, given by
\n$$
\nz\\in\\mathbb{C}\\mapsto\\frac{\\mathrm{e}^{-N|z|^2}}{\\pi}\\sum_{\\ell=0}^{N-1}\\frac{|\\sqrt{N}z|^{2\\ell}}{\\ell!}
\n$$
\nwhich is the analogue of the one given above for the Gaussian Unitary Ensemble. Using induction and integration by parts, it turns out that this density can be rewritten as
\n$$
\nz\\in\\mathbb{C}\\mapsto\\int_{N|z|^2}^\\infty\\frac{u^{N-1}\\mathrm{e}^{-u}}{\\pi(N-1)!}\\mathrm{d}u
\n=\\frac{\\Gamma(N,N|z|^2)}{\\pi}
\n$$
\nwhere $\\Gamma$ is the normalized incomplete Gamma function and where we used the identity
\n$$
\n\\mathrm{e}^{-r}\\sum_{\\ell=0}^n\\frac{r^\\ell}{\\ell!}=\\frac{1}{n!}\\int_r^\\infty u^n\\mathrm{e}^{-u}\\mathrm{d}u.
\n$$
\nNote that the function $t\\mapsto 1-\\Gamma(N,t)$ is the cumulative distribution function of the Gamma distribution with shape parameter $N$ and scale parameter $1$. As a consequence, denoting $X_1,\\ldots,X_N$ i.i.d. exponential random variables of unit mean, the density can be written as
\n$$
\nz\\in\\mathbb{C}\\mapsto\\frac{1}{\\pi}\\mathbb{P}\\left(\\frac{X_1+\\cdots+X_N}{N}\\geq|z|^2\\right).
\n$$
\nNow the law of large numbers implies that this density converges pointwise as $N\\to\\infty$ towards the density of the uniform distribution on the unit disk of the complex plane
\n$$
\nz\\in\\mathbb{C}\\mapsto\\frac{\\mathbf{1}_{|z|\\leq1}}{\\pi}.
\n$$
\nOne can also use i.i.d. Poisson random variables $Y_1,\\ldots,Y_N$ of mean $|z|^2$, giving
\n$$
\nz\\in\\mathbb{C}\\mapsto\\frac{1}{\\pi}\\mathbb{P}\\left(\\frac{Y_1+\\cdots+Y_N}{N}<1\\right),
\n$$
\nsee this former post<\/a>. Here are the plots with respect to $|z|$, followed by the Julia code.<\/p>\n
\"\"<\/p>\n
# Pkg.add(\"Distributions\") # https:\/\/juliastats.github.io\/Distributions.jl\r\n# Pkg.add(\"Plots\") # http:\/\/docs.juliaplots.org\/latest\/tutorial\/\r\nusing Plots\r\nusing Distributions \r\nN = 8\r\nX = collect(0:.01:3)\r\nY = ccdf(Gamma(N,1),N*X.^2)\/pi\r\nfunction unit(x)\r\n    if (x>1) return 0 else return 1\/pi end\r\nend # function\r\nY = [Y [unit(x) for x in X]]\r\npyplot()\r\nplot(X,Y,\r\n     title = @sprintf(\"Ginibre mean Empirical Spectral Distribution beta=2 N=%i\",N),\r\n     titlefont = font(\"Times\", 10),\r\n     label = [\"Ginibre mean ESD\" \"Equilibrium\"],\r\n     lw = 2,\r\n     linestyle = [:dash :solid],     \r\n     linecolor = [:darkblue :darkred],\r\n     grid = false,\r\n     border = false)\r\ngui()\r\n#savefig(\"ginibre.svg\")\r\nsavefig(\"ginibre.png\")<\/pre>\n
Further reading.<\/strong> The\u00a0Introducing Julia<\/a>\u00a0wikibook on Wikibooks<\/a>.<\/p>\n","protected":false},"excerpt":{"rendered":"
This short post is devoted to a couple of Julia programs. Gaussian Unitary Ensemble. It is the Boltzmann-Gibbs measure with density proportional to $$ (x_1,\\ldots,x_N)\\in\\mathbb{R}^N\\mapsto\\mathrm{e}^{-N\\sum_{i=1}^N x_i^2}\\prod_{i<j}(x_i-x_j)^2.…<\/p>\n
Continue readingPlaying a bit with Julia<\/span><\/a><\/div>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"iawp_total_views":198},"categories":[1],"tags":[],"_links":{"self":[{"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/posts\/10133"}],"collection":[{"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/comments?post=10133"}],"version-history":[{"count":42,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/posts\/10133\/revisions"}],"predecessor-version":[{"id":10217,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/posts\/10133\/revisions\/10217"}],"wp:attachment":[{"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/media?parent=10133"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/categories?post=10133"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/tags?post=10133"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}