Message Boards

GROUPS:

Staff Picks Mathematics Calculus Equation Solving Graphics and Visualization Wolfram Language Optimization Modeling Numerical Computation

Rauan Kaldybaev

Spherical harmonics interpolation method that minimizes an energy function

Rauan Kaldybaev

Posted 22 days ago

354 Views

1 Reply

4 Total Likes

Follow this post

Spherical harmonics interpolation method that minimizes an energy function
by Rauan Kaldybaev

Interpolation on a sphere is an essential problem in engineering, computer graphics, acoustics, and other fields, not to mention its mathematical interest. One approach to spherical interpolation is through spherical harmonics - these functions are smooth, numerically convenient, and easy to analyze. One issue, however, is that interpolation with spherical harmonics generally produces and overdetermined system, as spherical harmonics of degree

form a “package” of size

2l+1

and can only meaningfully be viewed together as a basis for the space of homogeneous polynomials of degree

that satisfy Laplace’s equations. This essay presents a method of interpolation that resolves this overdetermination by using the excess degrees of freedom to minimize a given energy function. A wisely chosen energy function makes high-degree modes costly, thus discouraging extreme variation, and is zero for the constant function. The resulting interpolation evaluates exactly to the reference values. The root-mean-squared error was found to be very accurately (

=0.984

) described by the empirical law

4/3

, where

is the number of data points. The interpolation respects spherical symmetry, so it does not introduce coordinate bias. To assess the effectiveness of min-energy interpolation, it was compared against the more traditional regularization-based interpolation on the real-world example of illumination of a sphere by lamps. With 110 data points, min-energy interpolation was found to be roughly 4 times more accurate than regularization-based interpolation. Then, to test the speed of convergence, the error of the two interpolations was plotted for different values of

; min-energy interpolation again performed better. The computation time of both interpolation methods is approximately the same. Min-energy interpolation is great for interpolation of expensive-to-compute mathematical functions. However, because it evaluates exactly to the reference values, it does not work well with noisy data. Min-energy interpolation basically corresponds to the limit of Tikhonov regularization with

λ->0

and a properly chosen regularization matrix.

Out[]=

Qualitative introduction. An overdetermined model

One problem with using spherical harmonics for interpolation is that they come in “packages”. There are

(L+1)

spherical harmonics with

l<=L

, meaning that the number of basis functions in such interpolation can only be a pure square - 1, 4, 9, 16, 25, and so on. If we have 11 points to interpolate, we have to stick with 16 spherical harmonics, which leaves 5 unused degrees of freedom. And this is what makes spherical harmonics interpolation somewhat nontrivial. One possible thing to do would be to “throw away” the extra harmonics - that is, just set some of the coefficients to zero. But that would introduce coordinate bias. The set of spherical harmonics with a particular value of

form a basis for the space of eigenfunctions of the Laplacian with eigenvalue

-l(l+1)

. If we forcefully remove some of the harmonics, we would be imposing an artificial constraint that contradicts the symmetry of the problem. And this, in turn, can dramatically decrease accuracy. As an example, imagine that we are asked to fit the linear model

ψ[x,y]=

to the data

ψ[1,3]=0

ψ[6,15]=26

. There are 3 degrees of freedom and only 2 constraints. One option would be to just set

, which leads to

=2.6

. However, this would mean choosing the

-axis as the preferred direction, and having preferred directions just doesn’t feel right. From a more practical perspective, the value

=5.2

is rather extreme: a more realistic estimate for the norm of

would be

26

(6-1)

(15-3)

. The only right direction for

is from

(1,3)

(6,15)

. And how exactly do we find “the right direction” in a more general case? For

points on a sphere, the set of coefficients isn’t even a vector (it’s more similar to a multi-index tensor). One very natural thing to do is to minimize a certain “energy function” under the constraint that the interpolation must evaluate to the given values at the reference points. Given that the energy function is chosen wisely, the resulting interpolation will not contain any coordinate bias. One can think of this interpolation as an elastic sphere that tries to minimize its potential energy. Regularization is not needed in this interpolation. Provided that the energy function is chosen well, high-frequency modes are penalized sufficiently to prevent extreme “spikiness”. Moreover, while regularization prevents an interpolation from evaluating exactly to the given values at the reference points, which invariably introduces inaccuracies, the minimization of an energy function happens under the constraint that the function has exactly the right values at the reference points. The lack of need for regularization also simplifies the interpolation process significantly, as the user does not have to choose the regularization parameters.

Mathematical description

Notation and conventions

The interpolation can be written as

ψ[ϕ,θ]=

(L+1)

∑

j=1

[ϕ,θ]

(

)

where

is an integer known as the truncation number;

-π<ϕ<=π

and

-π/2<=θ<=π/2

are the azimuthal and longitudinal angles, respectively;

are the spherical harmonics numbered by a single number such that

(the constant function) and the others can go in any order;

are the unknown coefficients. The formula (1) describes all possible linear combinations of the spherical harmonics

with

l<=L

; the upper limit of summation

(L+1)

comes from adding up the number of harmonics with different

’s:

1+3+5+...+(1+2L)=

(L+1)

. The requirement that the interpolation evaluates to the given values at the reference points can be written in matrix form as

Ax=b

(

)

where

is the (known) value of the interpolated function at the

-th point.

[

]

is the value of the

-th spherical harmonic at the

-th point. (If

were square, we would call it the design matrix). Because the basis functions are different,

is full-rank. Let us denote

the number of data points and

(L+1)

the number of spherical harmonics used for interpolation. Then, the matrix

mn

. To ensure that the model is not underdetermined,

n>=m

. The spherical harmonics in this essay are defined as follows:

≡

Cos[kϕ] k P l [Sin[θ]]	k>=0
Sin[\|k\|ϕ] \|k\| P l [Sin[θ]]	k<0

(

)

where

is a normalization factor and

is the associated Legendre polynomial. The spherical harmonics as defined by (3) are eigenfunctions of the Laplace operator with eigenvalue

-l(l+1)

and of the operator

∂

with eigenvalue

; they are also orthonormal. The reason to replace

ϕ



with sines and cosines is to make the basis functions real (since this essay will only consider interpolating real functions). The numbering used in the essay is

-1

-2

-1

, and so on; but again, the order doesn’t really change anything, the only convention that must be kept is that

Resolving the overdetermined system

When a model has more parameters than constraints, one possible solution is to use the extra degrees of freedom to minimize a certain energy function. One way to do that is to express the

-component vector

in terms of a smaller

-component vector

such that any vector

computed from

would automatically minimize the energy. This way, we would have

unknown components of

to solve exactly

equations. To make the minimization equations linear, let us only consider quadratic energy functions. If

ψ=

, then energy is given by

(

)

where

is a symmetric matrix of weights. Thus, the problem is minimizing the energy (4) while obeying the equation

Ax=b

. In other words, if the coefficient vector is

, then the value of energy must not change (to first order) if the coefficient vector were changed to

x->x+εδ

, with

ε->0

, for any vector

allowed by the constraint

Ax=b

. Mathematically, the gradient of

is orthogonal to every vector in the nullspace of

. Since the orthogonal complement of the nullspace is the rowspace, the gradient of

must lie in the rowspace of

. The same conclusion can be reached using Lagrange’s multipliers.

x∈R[A]

(

)

Given that most physical and real-world functions are defined up to a constant - for example, temperature - the energy function must not punish constant functions. In other words, adding

(a constant) to

should not change energy. What this means is that the first column and row of

should be zero.

then has the (somewhat problematic) form

=

0	0
0	W



(

)

The equation above is written in block matrix notation; here, the top-left zero is just the number 0, the off-diagonal zeros denote blocks of

n-1

zeros, and the matrix

is the “inner part” of

. Other than the constant function, there is no function that can be added to

such that

does not change. Correspondingly, the nullspace of

is just the zero vector, which means the matrix is invertible; this is going to be important. The equation (5) becomes



0	0
0	W

x∈R[A]

(

)

We see that the first entry

of the vector

does not play any role in the equation (7). A natural thing to do is to separate the vector into two components:

x=



(

)

where block matrix notation is used:

is a number and

is an

n-1

component vector. The equation (7) then becomes



∈R[A]

, or, equivalently,

(

)

for some

-component vector

. The matrix

was defined as

[

]

- in other words, the

-th column of

lists the values of the

-th spherical harmonic at the reference points.

, in particular, lists the values of

; but

is a constant, so

is the same for all

. A convenient way to write this is as follows:

A=(

)

(

)

where

is an

m1

column vector of the same number

and

is an

m(n-1)

matrix. The equation (9) then becomes

(

)

By looking at the first row of the equation, we get

v=0

. One way to achieve this equality is to declare

for

i>=2

and

∑

i=2

for some other vector

. In terms of the vector

, the second part of (11) becomes

Wy=

0	-1	-1	-1	...
	1
		1
			1
				...

(

)

where the unnamed matrix is of size

mm

and all the unfilled slots are zero; it can be written in piecewise notation as

[unnamedmatrix]

. The first component of

is not used here, so we can use it to set

in the formula (8). Thus, we have

x=Mu

(

)

where

M=

1	0
0	F



(

)

In the formula (14), the 1 in the top left corner is the number 1,

is an

(n-1)(m-1)

matrix, and the zeros represent blocks of zeros of appropriate size. The formula for

is as follows:

-1

-1	-1	-1	...
1
	1
		1
			...

(

)

By components,

(

-1

)

i,j+1

(

-1

)

i,1

. The unnamed matrix is of size

m(m-1)

; it is like the unnamed matrix in (12) but with the first column removed. It can be shown that whatever

is passed into the formula (13), the result is such that

x∈R[A]

. We have effectively taken care of the minimization part of the problem and now only need to satisfy

Ax=b

by finding the right

. The vector

has

components - this is just enough to satisfy the

constraints. Substituting

x=Mu

into

Ax=b

, we have

AMu=b

(

)

can be written as

(

)

, where

denotes an

n1

vector of

. Because

is symmetric,

-1

)

; denoting

C≡B

-1

, we can efficiently compute

j+1,i

1,i

. Finally, the vector

can be computed by substituting

into (13). It solves

Ax=b

and uses the remaining degrees of freedom to minimize

Numerical check


Relation to Tikhonov regularization

The formula (16) basically minimizes the quantity

(Ax-b)

(Ax-b)+λ

, where

λ->

. This can be seen from a physical metaphor of a small bead inside

Spherical harmonics interpolation method that minimizes an energy function

Qualitative introduction. An overdetermined model

Mathematical description

Notation and conventions

Resolving the overdetermined system

Numerical check

Relation to Tikhonov regularization

Numerical check
