Expected number of lines meeting four given lines or "what is 1.72..."

Given four random lines in $$\mathbb{R}P^3$$, how many lines intersect all of those lines? In the recent paper Probabilistic Schubert Calculus, Peter Bürgisser and Antonio Lerario discuss this question and much more general versions of it.

In Proposition 6.7, they determine the expected number of lines meeting four random lines in $$\mathbb{R}P^3$$ to be:

$$\operatorname{edeg}G(2,4) =\\ 2^{-13}\int_{[0,2\pi]^6}\left| \det\begin{pmatrix} \sin{t_1}\sin{s_1}&\sin{t_2}\sin{s_2}&\sin{t_3}\sin{s_3}\\ \cos{t_1}\sin{s_1}&\cos{t_2}\sin{s_2}&\cos{t_3}\sin{s_3}\\ \sin{t_1}\cos{s_1}&\sin{t_2}\cos{s_2}&\sin{t_3}\cos{s_3} \end{pmatrix}\right|dt_1dt_2dt_3ds_1ds_2ds_3$$

The integrand can be expanded to $$\left|\cos(s_2)\sin(s_1)\sin(s_3)\sin(t_2)\sin(t_1 - t_3)- \sin(s_2)\big(\cos(s_1)\sin(s_3)\sin(t_1)\sin(t_2 - t_3) + \cos(s_3) \sin(s_1)\sin(t_3)\sin(t_1 - t_2)\big)\right|$$

The absolute value in the integrand and the integration over six variables seem to make it difficult to evaluate this integral to high precision. I am interested in the exact value of this number. Bürgisser and Lerario computed it to be $$1.72$$ rounded to two digits; I ran $$10^{11}$$ Monte Carlo evaluations to obtain $$1.726225\dots$$ with an estimated error of $$7.3\cdot 10^{-6}$$

My questions are:

• Can the exact value of $$\operatorname{edeg}G(2,4)$$ be determined; is there perhaps a relation with other constants; is it, for example, algebraic over $$\mathbb{Q}[\pi]$$?
• If this is too difficult, what are ways of evaluating the above integral to higher accuracy numerically?

Update: In AdamP.Goucher's great answer, he provides more digits ($$1.726230876$$) and also uses a reformulation described by Matt F. in a comment to get an integrand without any absolute values:

$$G(2,4) = 2^{-6} \int \dfrac{h\left((y-u)^2+(v-x)^2\right)^2}{\left((1+(u + (y-u)c - (v-x)h)^2)(1+(x + (v-x)c + (y-u)h)^2)\prod_{\alpha\in\{u,y,x,v\}}(1+\alpha^2)\right)^{3/2}}$$ where we integrate over $$\mathbb{R}^5\times \mathbb{R}^+$$ (and $$h$$ is the positive variable).

Perhaps this form could come in handy for evaluating this integral to higher accuracy?!

Update 2 The new formula from L.Mathis and Antonio Lerario is very useful for calculating digits! The following mpmath code can returns in less than $$3$$ minutes

$$1.7262312489219034885256331685361697650475579915479447$$

(the last few digits might not be accurate) I expect to make that even much faster when solving the two integrals $$F$$ and $$G$$ symbolically first.

import functools
from mpmath import mp
@functools.lru_cache(maxsize=1000)
def F(u):
return mp.quad(lambda phi: (u*mp.sin(phi)**2)/(mp.sqrt(mp.cos(phi)**2 + u**2*mp.sin(phi)**2)), [0, mp.pi/2])
@functools.lru_cache(maxsize=1000)
def G(u):
return mp.quad(lambda phi: (mp.sin(phi)**2)/(mp.sqrt(mp.sin(phi)**2 + u**2*mp.cos(phi)**2)), [0, mp.pi/2])
def H(u):
return F(u)/G(u)
def L(u):
return F(u)*G(u)
def integrand3(u):
return L(u)**2*(1/H(u) - H(u))*mp.diff(H, u)/H(u)

dps = 50
mp.dps = dps

%time z = 3*mp.quad(integrand, [0,1]); z

• Note that the integrand is periodic in each direction with period pi, so you can cut your amount of work by a factor of 64. There's another factor of 6 from the symmetry of permuting (s1, t1), (s2, t2) and (s3, t3). A preliminary run using Gauss-Legendre quadrature with ~ 99^6/6 points gave an approximation of 1.726209603407 within 12 seconds; I'll try a larger net (maybe 199^6) reasonably soon. Jan 26 '17 at 19:47
• 1.726228333055 using 199^6 points. Now trying 299^6 points... Jan 26 '17 at 20:07
• @AdamP.Goucher: I didn't use a grid, but choose the points uniformly random, so I guess it doesn't make much of a difference?! Jan 26 '17 at 21:51
• I'm using Gauss-Legendre quadrature, so the grid is non-uniform. It's unlikely (but not impossible) to be any worse than Monte Carlo, and can sometimes be significantly better (e.g. for smooth functions). Also, this works out faster since I don't need to generate random numbers or perform trig functions on the fly. By the way, my latest result is 1.726230356553 with 299^6 points (in 95 minutes). I'm doing an overnight run with 401^6 points, starting now... Jan 26 '17 at 21:58
• @AdamP.Goucher: great, thanks! that should give a few more digits.. Jan 26 '17 at 22:03

The integrand is periodic modulo $\pi$ in each variable, so it suffices to integrate each variable over $[0, \pi]$ and replace the constant factor by $2^{-7}$.

If we were to apply a change of variables (e.g. set $x = \cos(s_1)$ and similarly for the other five variables), we would have an integral of a piecewise-algebraic function which thus belongs to Kontsevich and Zagier's 'ring of periods'. Alas, there is no proven algorithm for determining whether periods are expressible in terms of elementary functions.

Now, using Gauss-Legendre integration with $N$ points in each variable, we can approximate it with a weighted sum of $N^6$ evaluations of the integrand. Since the absolute value of the determinant is unchanged under permuting the columns, and zero whenever two columns are equal, we can reduce this to $\binom{N^2}{3}$ function evaluations.

I managed to implement this in C code, where each iteration, amortized, only takes four floating-point additions, three multiplications, and an absolute-value calculation. With several hours running on a 144-core machine, it could do the calculation for both $N = 401$ and $N = 409$. The results were 1.726230867 and 1.726230885, respectively.

Consequently, the integral is roughly 1.726230876, with an expected error on the order of $10^{-8}$.

EDIT: Matt F. posted an algebraic form for the integral, by taking $u = \cot(t1)$, $v = \cot(s2)$, $w = \cot(t3)$, $x = \cot(s1)$, $y = \cot(t2)$, $z = \cot(s3)$:

$$I = 2^{-7} \int \dfrac{|uv-vw+wx-xy+yz-zu|}{[(1+u^2)(1+v^2)(1+w^2)(1+x^2)(1+y^2)(1+z^2)]^{3/2}}$$

where the integral is taken over $(u,v,w,x,y,z) \in \mathbb{R}^6$. Observing that the numerator is twice the area of a triangle with vertices $(u, x), (y, v), (w, z)$, we can perform a change of variables by setting:

$$w = u + (y-u)c - (v-x)h \textrm{ ; } z = x + (v-x)c + (y-u)h$$

where $c \in \mathbb{R}$ and $h \in \mathbb{R}^{+}$. Taking the Jacobian into account, the integral becomes:

$$I = 2^{-6} \int \dfrac{h[(y-u)^2+(v-x)^2]^2}{[(1+u^2)(1+y^2)(1+(u + (y-u)c - (v-x)h)^2)(1+x^2)(1+v^2)(1+(x + (v-x)c + (y-u)h)^2)]^{3/2}}$$

where $(u,y,x,v,c,h) \in \mathbb{R}^5 \times \mathbb{R}^+$. Although the integral is uglier, it satisfyingly has no absolute value operator.

• A perspicuous algebraic form for the integral is: $$2^{-13}\int_{[-\infty,\infty]^6} \frac{\big|uv-vw+wx-xy+yz-zu \big|}{\big( (1+u^2)(1+v^2)(1+w^2)(1+x^2)(1+y^2)(1+z^2) \big)^{3/2}}du\,dv\,dw\,dx\,dy\,dz$$ which uses $u=\cot(t1),\ v=\cot(s2),\ w=\cot(t3),\ x=\cot(s1),\ y=\cot(t2),\ z=\cot(s3)$ Jan 29 '17 at 17:24
• nice work, Adam! @MattF.: How did you find this algebraic form? Jan 30 '17 at 12:03
• @MoritzFirsching: first factor out two sins from each column, leaving only 1's and cot's in the determinant; then write sin as 1/sqrt(1+cot^2); and then multiply by some more of those expressions to account for dt1/du, etc. Finally make the coefficient 2^-7, since u runs through the real line twice as t1 goes around the circle. Jan 30 '17 at 17:54
• @MattF. Presumably the $2^{-13}$ is meant to be $2^{-7}$, then? Jan 31 '17 at 0:23
• @AdamP.Goucher That would be great! Jan 31 '17 at 9:10

If anyone is still interested, Antonio Lerario and I recently published a paper: Probabilistic Schubert Calculus: Asymptotics, in which we give a more convenient formula to compute this number.

What was called before $$\operatorname{edeg}G(2,4)$$ is now denoted by $$\delta_{1,3}$$ and we give a line integral formula in Proposition 24:

\begin{aligned} \delta _{1,3}=-6 \pi ^{4}\int _{0}^1 L(u)^{2}\mathrm {sinh}(w(u))w'(u)\mathrm{d}u \end{aligned} where $$L=F⋅G$$ and $$w=log(F/G)$$ with \begin{aligned} F(u)&:=\int _{0}^{\pi /2} \frac{u \ \sin ^2(\theta )}{\sqrt{\cos ^2\theta +u^2 \sin ^2\theta }}\mathrm {d}\theta \\ G (u)&:=\int _{0}^{\pi /2} \frac{\sin ^2(\theta )}{\sqrt{\sin ^2\theta +u^2 \cos ^2\theta }}\mathrm {d}\theta . \end{aligned}

We didn't run any advanced numerical evaluation but this looks much nicer than the previous formula. I hope this will help. It is still unknown to us if this number can be expressed as a closed formula using special functions.

• Bravo! I'm still interested in seeing if this is actually easier to evaluate numerically.. Nov 4 '20 at 10:21
• In the linked paper, in Proposition 24, there is no "$u$" in the numerator of the integrand in the definition of $G(u)$. Is this a typo? Nov 4 '20 at 10:52
• I'm trying to use your formula, so far without much luck gist.github.com/mo271/649b781bb362931de9483c41214345eb I can calculate the integral, but I get 9.173... Nov 4 '20 at 11:29
• Maple evaluates $F(u)$ and $G(u)$ in terms of elliptic functions. Nov 4 '20 at 11:33
• Thank you for pointing it out! The correct form should be in Proposition 24 (without the "$u$" in the numerator of $G$). I edited the answer. I hope this will give better results. Nov 5 '20 at 11:03