High Order Numerical Cubature (Clean and Simple)
Forget every article I sent you about numerical quadrature/cubature ...
last night, in the night, vast simplification came to me ... it is much
easier to generalize the midpoint rule than the trapezoid rule for cubature
in arbitrarily many dimensions.
There is historical precedence for this idea ... ... the
*Gragg–Bulirsch–Stoer
(GBS) * numerical solution method for ODE's uses a combination of the
midpoint rule and Richardson extrapolation for high order accuracy. We are
going to generalize this to multivariate quadrature or "cubature."
Given a simplex (the convex hull of a finite set of points having affine
independence) we form a tree (more of a dendrite, really) whose root node
is the average of the vertices of the simplex (i.e. the "barycenter" of the
simplex), and whose descendants become denser in the interior of the
simplex with each successive generation.
This "space filling" has to be done with sufficient fractal self
symmetrical regularity to ensure that our method can be highly accurate and
efficient. Don't worry ... it turned out to be way simpler than I could
have possibly hoped for!
As we said, the first node of the tree is the baryceenter of the simplex.
The immediate children of that node are the points precisely half way from
the baryceenter to the respective vertices of the simplex. [The simplex is
the convex hull of these vertices.]
The respective arrows (or "vectors") from the barycenter to the vertices
serve as the only directions needed for finding all of the children of
every node.
The only question is "how far?"
Answer: the distance from a parent to a child is half as far as it was in
the previous generation ... i.e. half as far as from grandparent to parent.
Now, let g(k) represent the set of nodes in the kth generation down from
the root node. Then g(1) and g(2), respectively, represent the children and
grandchildren of the barycenter, while g(0) consists of only the barycenter
itself.
Also let G(k) be the union of all of the generations from g(0) to g(k),
inclusive.
Now let f be a real valued function defined at every node of our tree.
And let A(k) be the average of all of the f values of all of the nodes of
the set G(k), i.e. from the barycenter value down to all of the values of
the nodes of the k-th generation.
If f is a reasonable function, then as k increases, A(k) approaches the
integral of f over the simplex, divided by its "capacity" (length, area,
volume, etc ... depending on the dimension of the simplex).
Let's denote this limiting value of A(k) by the Greek letter Alpha.
In particular, when f is a density function, Alpha is the average density
of the simplex, because the integral of density is mass, and mass divided
by capacity is average density ... in units of g/(cm)^d or kg/m^d, where d
is the dimension of the simplex.
Similarly, if f represents temperature or pressure, Alpha will be the
average temperature or pressure experienced by the simplex.
If you need to know the value of the integral of f over the simplex, just
multiply its capacity by Alpha.
To get the capacity of a simplex, let M be a matrix whose rows are the
respective displacement vectors from one fixed vertex v0 to the others ...
v1-v0, v2-v0, etc.
The capacity of a d-dimensional simplex is the square root of the
determinant of the matrix product of M and its transpose, divided by
d-factorial :
Sqrt(det MM')/d!,
where M' is the transpose of M.
That's everything ... except the Richardson extrapolation to increase the
order of convergence from second order to order 2n, where n depends on the
number of Richardson transformations of the sequence [see penultimate
paragraph below]:
We start with the sequence
A(0), A(1), .... A(j), ...
Then transform the A sequence to a B sequence, where
B(0)=(4A(1)-A(0))/3 ...
B(j)=(4A(j+1)-A(j))/3 ...
Then transform the B's to a sequence of C's:
C(j)=(16B(j+1)-B(j))/15 ...
.......
D(j)=(4^3*C(j+1)-C(j))/(4^3-1)
etc.
The suggested order is ..
A0), A(1),B(0), A(2), B(1), C(0), A(3),B(2), C(1), D(0), etc....
At each step go as far down the alphabet as possible, given the previous
steps.
When you reach a new letter, i.e. the zero_th term in a new sequence,
compare it with the most recent value of the preceding letter to see if the
value is still progressing, or if it has reached a point of diminishing
returns.
Test it on polynomials: by the time the extrapolation reaches the n_th
letter of the alphabet, any polynomial function f of degree 2n or less
should be averaged exactly (up to roundoff error).
If your polynomials have integer coefficients, and the coordinates of the
vertices are integers or precise fractions, then there is no excuse for
roundoff error ... the answers must have infinite precision! (unless you
cannot afford it)
Have Fun!
Forest
On Thu, Feb 3, 2022 at 12:46 PM Forest Simmons forest.simmons21@gmail.com
wrote:
Forget every article I sent you about numerical quadrature/cubature ...
last night, in the night, vast simplification came to me ... it is much
easier to generalize the midpoint rule than the trapezoid rule for cubature
in arbitrarily many dimensions.
There is historical precedence for this idea ... ... the *Gragg–Bulirsch–Stoer
(GBS) * numerical solution method for ODE's uses a combination of the
midpoint rule and Richardson extrapolation for high order accuracy. We are
going to generalize this to multivariate quadrature or "cubature."
Given a simplex (the convex hull of a finite set of points having affine
independence) we form a tree (more of a dendrite, really) whose root node
is the average of the vertices of the simplex (i.e. the "barycenter" of the
simplex), and whose descendants become denser in the interior of the
simplex with each successive generation.
This "space filling" has to be done with sufficient fractal self
symmetrical regularity to ensure that our method can be highly accurate and
efficient. Don't worry ... it turned out to be way simpler than I could
have possibly hoped for!
As we said, the first node of the tree is the baryceenter of the simplex.
The immediate children of that node are the points precisely half way from
the baryceenter to the respective vertices of the simplex. [The simplex is
the convex hull of these vertices.]
The respective arrows (or "vectors") from the barycenter to the vertices
serve as the only directions needed for finding all of the children of
every node.
I think you are onto something. The exact way you phrased the algorithm
seems to only create points along the lines joining a vertex and the
barycenter. If I understood that correctly, then it is not "space filling".
But we can easily resolve that. Your key insight is that the simplex is
convex and a sum of vertices and/or points interior to the simplex will
give you a point that is also interior to the simplex.
One way to create points that are guaranteed to be inside the simplex is to
grab the vertices and do any weighted sum, as long as the weights add up to
- You could just pick those points randomly, but that wouldn't generalize
the midpoint rule the way that the barycenter does. Finally, everything
that we've said up to this point applies to any convex hull; not just
simplices. In an N-dimensional space you can grab any set of N+1 vertices
from your convex hull and their barycenter will be somewhere inside the
convex hull as well.
This last comment about taking N+1 vertices does not fundamentally resolve
anything. That's still a small number of points. To truly get a series of
space-filling points that generalizes the midpoint rule, we can keep
computing barycenters:
- Let S = some simplex. S has N+1 points in an N-dimensional space.
- Compute b = Barycenter(S)
- Add 'b' to the list of points: S -> {S, b}
- Select a random set of N+1 points from S
- Compute their barycenter and add that to S.
- GOTO (4)
The same algorithm can work for a convex hull.
- Let C = some convex hull with > N+1 points in an N-dimensional space.
- Select a random set of N+1 points from C.
- Compute the barycenter of those points and add that to C.
- GOTO (2)
Does this sound right or did I just misunderstand your algorithm?
Cheers,
Dr. Daniel Carrera
Postdoctoral Research Associate
Iowa State University
This is a slight revision of the original with some clarifications to
prevent all posible misunderstandings! (Sure)
In particular, if G is a set of points and V is a set of vectors, let G+V
be the set obtained by displacing each point of G by every vector of V:
G+V={g+v| (g,v) in G×V)},
where G×V is the Cartesian product of G and V.
So in general (when there are no double representations of points) the
number of points in this kind of sum will be the same as #(G×V), wich is
just (#G)*(#V).
Given a simplex (the convex hull of a finite set of points having affine
independence) we form a tree whose root node is the average of the vertices
of the simplex (i.e. the "barycenter" of the simplex), and whose
descendants become denser in the interior of the simplex with each
successive generation.
This "space filling" has to be done with sufficient fractal self
symmetrical regularity to ensure that our method can be highly accurate and
efficient. Don't worry ... it turned out to be way simpler than I could
have possibly hoped for!
As we said, the first node of the tree is the baryceenter of the simplex.
The immediate children of that node are the points precisely half way from
the barycenter to the respective vertices of the simplex. [The simplex
itself is the convex hull of these vertices.]
The respective vectors from the barycenter to the vertices form a set V of
displacement vectors that serve as the only directions needed for finding
all of the children of every node.
Once we know the directions, the only question that remains is "how far?"
(in those directions).
Answer: the distance from a parent to a child is half as far as it was in
the previous generation ... i.e. half as far as from grandparent to parent.
Now, let G(k) represent the set of nodes in the k_th generation down from
the root node. Then G(1) and G(2), respectively, represent the children and
grandchildren of the barycenter, while G(0) consists of only the barycenter
itself.
In summary, for each natural number k, the k_th generation is obtained from
the previous generation by the recursive formula G(k)=G(k-1)+V/2^k .
Now let f be a real valued function defined at every node of our tree.
And let A(k) be the average of all of the f values of all of the nodes of
the set G(k), i.e.
A(k)=Sum{f(g)| g in G(k)}/#G(k).
If f is a reasonable function, then as k increases, A(k) approaches the
average value of f on the simplex, which is defined by calculus as the
integral of f over the simplex, divided by the d-dimensional capacity of
the simplex. The "capacity" of a simplex is its length, area, volume, etc
... depending on the dimension d of the simplex.
Let's denote the limiting value of A(k) by the Greek letter Alpha(f).
In particular, when f is a density function, Alpha(f) is the average
density of the simplex, because the integral of density is mass, and mass
divided by capacity is average density ... in units of g/(cm)^d or kg/m^d,
where d is the dimension of the simplex.
Similarly, if f represents temperature or pressure, Alpha(f) will be the
average temperature or pressure experienced by the simplex.
If you need to know the mass of the simplex (when f is its density
function), just multply its average density Alpha(f) by its capacity.
To get the capacity of a simplex, let M be a matrix whose rows are the
respective displacement vectors from one fixed vertex v0 to the others ...
v1-v0, v2-v0, etc.
The capacity of the d-dimensional simplex is the square root of the
determinant of the matrix product of M and its transpose M', divided by d.
In other words...
Capacity=Sqrt(det MM')/d
When M is a square matrix this capacity computation is simply
Capacity=|det M|/d.
That's everything ... except the Richardson extrapolation to increase the
order of convergence from second order to order 2n, where n depends on the
number of Richardson transformations of the sequence [for more information
on the rate of convergence or order of accuracy, see the last remark below]:
We start with the sequence
A(0), A(1), .... A(j), ...
Then transform the A sequence to a B sequence, where ...
B(0)=(4A(1)-A(0))/3 ...
B(j)=(4A(j+1)-A(j))/3 ...
Then transform the B's to a sequence of C's:
C(j)=(16B(j+1)-B(j))/15 ...
.......
D(j)=(4^3*C(j+1)-C(j))/(4^3-1)
etc.
The suggested order is ..
A0), A(1),B(0), A(2), B(1), C(0), A(3),B(2), C(1), D(0), etc....
At each step go as far down the alphabet as possible, given the results
needed from the previous steps.
When you reach a new letter, i.e. the zero_th term in a new sequence,
compare it with the most recent value of the preceding letter to see if the
approximation is still progressing, or if it has reached a point of
diminishing returns.
Test it on polynomials: by the time the extrapolation reaches the n_th
letter of the alphabet, any polynomial function f of degree 2n or less
should be averaged quite precisely (up to roundoff error).
-Forest
Hello Forest,
On Sat, Feb 5, 2022 at 2:00 PM Forest Simmons forest.simmons21@gmail.com
wrote:
This is a slight revision of the original with some clarifications to
prevent all posible misunderstandings! (Sure)
In particular, if G is a set of points and V is a set of vectors, let G+V
be the set obtained by displacing each point of G by every vector of V:
G+V={g+v| (g,v) in G×V)},
where G×V is the Cartesian product of G and V.
So in general (when there are no double representations of points) the
number of points in this kind of sum will be the same as #(G×V), wich is
just (#G)*(#V).
Aha! I think this added the clarity that I needed. Now I understand what
you're saying. To make sure we are on the same wavelength, I wrote a script
to plot a schematic of how I understand your algorithm to work, followed by
the first few generations up to G(6). Here is the plot:
Unfortunately, it doesn't seem to be space-filling. It looks like you have
reinvented Sierpiński's gasket.
Let me know if I misunderstood the algorithm again.
Cheers,
Dr. Daniel Carrera
Postdoctoral Research Associate
Iowa State University
Nice graphics! We should call it "shingles".
I wonder if you could do the same thing for an equilateral triangle ...
say with vertices at (-1,0), (1, 0) and (0, sqrt(3)).
El sáb., 5 de feb. de 2022 1:55 p. m., Daniel Carrera dcarrera@gmail.com
escribió:
Hello Forest,
On Sat, Feb 5, 2022 at 2:00 PM Forest Simmons forest.simmons21@gmail.com
wrote:
This is a slight revision of the original with some clarifications to
prevent all posible misunderstandings! (Sure)
In particular, if G is a set of points and V is a set of vectors, let G+V
be the set obtained by displacing each point of G by every vector of V:
G+V={g+v| (g,v) in G×V)},
where G×V is the Cartesian product of G and V.
So in general (when there are no double representations of points) the
number of points in this kind of sum will be the same as #(G×V), wich is
just (#G)*(#V).
Aha! I think this added the clarity that I needed. Now I understand what
you're saying. To make sure we are on the same wavelength, I wrote a script
to plot a schematic of how I understand your algorithm to work, followed by
the first few generations up to G(6). Here is the plot:
Unfortunately, it doesn't seem to be space-filling. It looks like you have
reinvented Sierpiński's gasket.
Let me know if I misunderstood the algorithm again.
Cheers,
Dr. Daniel Carrera
Postdoctoral Research Associate
Iowa State University
On Sat, Feb 5, 2022 at 3:55 PM Daniel Carrera dcarrera@gmail.com wrote:
Unfortunately, it doesn't seem to be space-filling. It looks like you have
reinvented Sierpiński's gasket.
Incidentally, my own idea of using multiple barycenters doesn't work
either. It gives points very clustered in the center of the simplex.
I did find a solution that works great for a 2D simplex, but I don't see an
obvious way to generalize it. We begin by generating a pair of points (u,v)
from the 2D Sobol sequence. We want to map these to (x,y) inside the
simplex in a way that retains at least some of the uniformity guarantees of
the Sobol sequence. The simplex is defined by the boundary:
x + y < 1
That means for x = 0 there are more possible values of y than for x = 0.9,
so we should produce more values close to x=0. This is a familiar problem
in the context of drawing random numbers and it has a known solution.
Let y = f(x) = 1 - x be (not normalized) the probability of generating
the value x. The cumulative distribution is:
F(x) = x - x^2/2
Normalize it to F(x) = 2x - x^2 so that we have the property that F(0) = 0 and F(1) = 1. So we can let F(x) be equal to the first Sobol number:
F(x) = u
=> x = 1 - sqrt(1 - u)
This gives an x with the correct probability distribution and we are free
to draw y randomly between 0 and (1-x). We accomplish that with the
second Sobol number:
=> y = (1 - x)*v
The result is a nicely uniform distribution of points without the
clustering that you'd expect from random points:
There's probably a way to generalize these ideas to higher dimensions.
Cheers,
Dr. Daniel Carrera
Postdoctoral Research Associate
Iowa State University
Ok. I think I have a general solution of the standard n-dimensional simplex.
Let {x1, x2, ... , xn} be the coordinates.
For any given value of x1 the remaining coordinates can go from 0 to (1-x1)
and span an (n-1)-simplex with volume V ~ (1-x1)^(n-1). To distribute the
points uniformly, the probability of choosing a value x1 should be
proportional to that volume:
Prob: f1(x1) ~ (1-x1)^(n-1)
The cumulative distribution is:
F1(x1) = Integral from 0 to x1 of f1(x1)
If we integrate this and normalize the integral so that F1(1) = 1 we get
F1(x1) = 1 - (1 - x1)^n
Let us now draw n values from the Sobol sequence. The way the Sobol
sequence works, you specify the number of dimensions and Sobol_nD returns a
single n-point from the unit N-cube.
{s1, s2, ..., sn} <-- Sobol_nD()
Since the cumulative distribution F1() goes from 0 to 1, I can set F1 = s1
or F1 = 1-s1. As it turns out, the latter gives a slightly nicer formula.
Let F1(x1) = 1 - s1
=> x1 = 1 - s1^(1/n)
Now that I have generated x1, I can generate other values recursively.
Suppose that I have generated values {x1, x2, ..., x{k-1}} and now I want
to generate xk. We have that 0 < xk < 1 - x{k-1}. For any choice of xk the
remaining coordinates span an (N-k)-simplex with volume (1-xk)^(N-k) so the
probability of choosing xk has to be proportional to that. So we setup the
probability distribution, integrate to get the cumulative distribution
Fk(xk), and normalize so that Fk(xk = 1 - x{k-1}) = 1. When I got to this
point the formula was slightly less ugly if I chose Fk(xk) = sk. Solve for
xk and I obtain a formula for xk based on the kth Sobol number and the
previous point x{k-1}.
The formula is slightly less ugly if I write it as a formula for x{k+1}
instead of xk. Long story short:
x1 = 1 - s1^(1/n)
x{k+1} = 1 - [ 1 - s{k+1} * ( 1 - xk^(n-k) ) ]^(1 / (n-k) )
As an example, let's compute the last coordinate xn. In this instance there
is a lot of cancellation and we get:
xn = sn * (1 - x{n-1})
I haven't tested this, but this should produce uniformly distributed points
{x1,x2,...,xn} inside the standard n-simplex.
Cheers,
Daniel
On Sat, Feb 5, 2022 at 10:32 PM Daniel Carrera dcarrera@gmail.com wrote:
On Sat, Feb 5, 2022 at 3:55 PM Daniel Carrera dcarrera@gmail.com wrote:
Unfortunately, it doesn't seem to be space-filling. It looks like you
have reinvented Sierpiński's gasket.
Incidentally, my own idea of using multiple barycenters doesn't work
either. It gives points very clustered in the center of the simplex.
I did find a solution that works great for a 2D simplex, but I don't see
an obvious way to generalize it. We begin by generating a pair of points
(u,v) from the 2D Sobol sequence. We want to map these to (x,y) inside the
simplex in a way that retains at least some of the uniformity guarantees of
the Sobol sequence. The simplex is defined by the boundary:
x + y < 1
That means for x = 0 there are more possible values of y than for x = 0.9,
so we should produce more values close to x=0. This is a familiar problem
in the context of drawing random numbers and it has a known solution.
Let y = f(x) = 1 - x be (not normalized) the probability of generating
the value x. The cumulative distribution is:
F(x) = x - x^2/2
Normalize it to F(x) = 2x - x^2 so that we have the property that F(0) = 0 and F(1) = 1. So we can let F(x) be equal to the first Sobol number:
F(x) = u
=> x = 1 - sqrt(1 - u)
This gives an x with the correct probability distribution and we are
free to draw y randomly between 0 and (1-x). We accomplish that with the
second Sobol number:
=> y = (1 - x)*v
The result is a nicely uniform distribution of points without the
clustering that you'd expect from random points:
There's probably a way to generalize these ideas to higher dimensions.
Cheers,
Dr. Daniel Carrera
Postdoctoral Research Associate
Iowa State University
--
Dr. Daniel Carrera
Postdoctoral Research Associate
Iowa State University
Good work!
That's much better than what I was thinking!
El dom., 6 de feb. de 2022 1:22 a. m., Daniel Carrera dcarrera@gmail.com
escribió:
Ok. I think I have a general solution of the standard n-dimensional
simplex.
Let {x1, x2, ... , xn} be the coordinates.
For any given value of x1 the remaining coordinates can go from 0 to
(1-x1) and span an (n-1)-simplex with volume V ~ (1-x1)^(n-1). To
distribute the points uniformly, the probability of choosing a value x1
should be proportional to that volume:
Prob: f1(x1) ~ (1-x1)^(n-1)
The cumulative distribution is:
F1(x1) = Integral from 0 to x1 of f1(x1)
If we integrate this and normalize the integral so that F1(1) = 1 we get
F1(x1) = 1 - (1 - x1)^n
Let us now draw n values from the Sobol sequence. The way the Sobol
sequence works, you specify the number of dimensions and Sobol_nD returns a
single n-point from the unit N-cube.
{s1, s2, ..., sn} <-- Sobol_nD()
Since the cumulative distribution F1() goes from 0 to 1, I can set F1 = s1
or F1 = 1-s1. As it turns out, the latter gives a slightly nicer formula.
Let F1(x1) = 1 - s1
=> x1 = 1 - s1^(1/n)
Now that I have generated x1, I can generate other values recursively.
Suppose that I have generated values {x1, x2, ..., x{k-1}} and now I want
to generate xk. We have that 0 < xk < 1 - x{k-1}. For any choice of xk the
remaining coordinates span an (N-k)-simplex with volume (1-xk)^(N-k) so the
probability of choosing xk has to be proportional to that. So we setup the
probability distribution, integrate to get the cumulative distribution
Fk(xk), and normalize so that Fk(xk = 1 - x{k-1}) = 1. When I got to this
point the formula was slightly less ugly if I chose Fk(xk) = sk. Solve for
xk and I obtain a formula for xk based on the kth Sobol number and the
previous point x{k-1}.
The formula is slightly less ugly if I write it as a formula for x{k+1}
instead of xk. Long story short:
x1 = 1 - s1^(1/n)
x{k+1} = 1 - [ 1 - s{k+1} * ( 1 - xk^(n-k) ) ]^(1 / (n-k) )
As an example, let's compute the last coordinate xn. In this instance
there is a lot of cancellation and we get:
xn = sn * (1 - x{n-1})
I haven't tested this, but this should produce uniformly distributed
points {x1,x2,...,xn} inside the standard n-simplex.
Cheers,
Daniel
On Sat, Feb 5, 2022 at 10:32 PM Daniel Carrera dcarrera@gmail.com wrote:
On Sat, Feb 5, 2022 at 3:55 PM Daniel Carrera dcarrera@gmail.com wrote:
Unfortunately, it doesn't seem to be space-filling. It looks like you
have reinvented Sierpiński's gasket.
Incidentally, my own idea of using multiple barycenters doesn't work
either. It gives points very clustered in the center of the simplex.
I did find a solution that works great for a 2D simplex, but I don't see
an obvious way to generalize it. We begin by generating a pair of points
(u,v) from the 2D Sobol sequence. We want to map these to (x,y) inside the
simplex in a way that retains at least some of the uniformity guarantees of
the Sobol sequence. The simplex is defined by the boundary:
x + y < 1
That means for x = 0 there are more possible values of y than for x =
0.9, so we should produce more values close to x=0. This is a familiar
problem in the context of drawing random numbers and it has a known
solution.
Let y = f(x) = 1 - x be (not normalized) the probability of generating
the value x. The cumulative distribution is:
F(x) = x - x^2/2
Normalize it to F(x) = 2x - x^2 so that we have the property that F(0) = 0 and F(1) = 1. So we can let F(x) be equal to the first Sobol number:
F(x) = u
=> x = 1 - sqrt(1 - u)
This gives an x with the correct probability distribution and we are
free to draw y randomly between 0 and (1-x). We accomplish that with the
second Sobol number:
=> y = (1 - x)*v
The result is a nicely uniform distribution of points without the
clustering that you'd expect from random points:
There's probably a way to generalize these ideas to higher dimensions.
Cheers,
Dr. Daniel Carrera
Postdoctoral Research Associate
Iowa State University
--
Dr. Daniel Carrera
Postdoctoral Research Associate
Iowa State University