\documentstyle[12pt]{article}
\begin{document}


{\bf Assessing and Improving Helioseismological Inferences of Solar Structure}\\

P.B. Stark, Principal Investigator\\
Department of Statistics and Space Sciences Laboratory\\
University of California\\
Berkeley CA 94720-3860\\


{\bf Abstract:}\\

Acoustic vibrations of the Sun contain information about the distributions of
soundspeed, material density, and angular velocity within the Sun.
Current methods for imaging these properties of the solar interior 
from helioseismic eigenfrequencies produce images that might resemble the Sun,
but the formal uncertainties of the techniques do not accurately quantify
possible differences between the models and the truth.
Nor can the methods incorporate physical constraints that might reduce the
uncertainty.
The effects of data errors, modeling algorithms, and theoretical 
approximations on the reliability of the images are largely unexplored; 
published uncertainties approximate the variability of the
estimation method conditional on a particular parametric representation of the
solar property, and a host of other assumptions, without regard for possible 
biases.
The proposed research is to develop and apply methods to study both the
variability and the bias of estimates of the distribution of soundspeed
in the solar interior and the distribution of angular velocity in the solar
interior, under less restrictive assumptions, incorporating prior physical
constraints such as bounds on energy or
rotationally-induced oblateness to control the bias.

Shorter, but still conservative,
confidence intervals for solar properties can be found
by using unusual functionals for measuring fit to the data, rather than
conventional $l_p$ misfit norms such as the common sum of squared
residuals. Such alternatives will be explored and applied to
find simultaneous confidence envelopes for the distribution of soundspeed in the
solar interior, incorporating uncertainties due to surface effects and to
neglected higher-order terms in the asymptotics, and to make direct inferences
about spatial variations in the distribution of angular velocity within the Sun
from the splitting of normal mode frequencies. 
A fringe benefit of the new misfit measures will be ``teleological'' data
compression techniques that enhance the ability to
image the solar interior from normal mode data.

Proposed Research.

1) Testing Hypotheses About Rotation.

Eigenmodes of the Sun with the same radial order and angular degree would have
the same eigenfrequency if the Sun were spherically symmetric.
Rotation of the Sun is one origin of asymmetry that ``splits''
degenerate eigenfrequencies.  To first order, the odd functional part
of the splitting of individual singlets from the central multiplet
frequency is a linear function of the angular velocity; asymmetry in the
acoustic velocity structure contributes to the even part of the
splitting. Observations of the splitting of normal mode frequencies can and have
been used to estimate the distribution of angular velocity within the
Sun.  It is common practice in these studies to take the soundspeed
distribution within the Sun as known; within that approximation, the
splitting is linearly related to the angular-velocity distribution, to first
order in the angular velocity.  Let
$\delta_j$ denote the $j$th observed frequency splitting, $\epsilon_j$ the
observational error in its measurement, $\Omega$ the distribution of angular
velocity within the Sun, and $B$ the Sun.  Within the customary approximation,
\begin{equation}
\delta_j = \int_B K_j({\bf r}) \Omega({\bf r}) d {\bf r} + \epsilon_j, \;\;j=1, \ldots,
n .
\end{equation}
It is often assumed (and reported by observers) that the errors 
$\epsilon_j$ are independent, zero mean
random variables with known variances $\sigma_j^2$.

2) Tailored Kernels by Numerical Optimization

Fitted models of solar rotation have many intriguing features.  However, since
the error analysis is not complete, it is not possible to infer from a feature
of a fitted model whether or not the Sun is likely to possess the same feature.
Genovese, Stark and Thompson (Ap.J., 1995) show that from Big Bear Solar
Observatory (BBSO) data (Woodard and Libbrecht (Ap.J., 1993), even
with very restrictive assumptions, inferences based on pointwise estimates of 
the rotation rate can not even rule out solid-<body> rotation. To make sound
inferences requires abandoning the search for entire rotation models.

Gough, Sekii and Stark (Ap.J., 1995, submitted) present a variety of ways to
make simultaneous inferences about a collection of linear functionals of solar
structure, and illustrate a promising approach I propose to pursue:
to make inferences about the solar interior by constructing functionals that
are directly sensitive to hypotheses of interest.  One may then test the
hypotheses by finding confidence intervals for the functionals.

For example,
suppose we wish to determine whether or not the average rotation rate increases
with depth across the convection zone.  
For a given set of $n$ real numbers $\{\beta_j \}_{j=1}^n$, define the linear
combination
\[
f_\beta({\bf r}) \equiv \sum_{j=1}^n \beta_j K_j ({\bf r} ) .
\]
Let $B^-$ and $B^+$ be non-overlapping spherical shells such that the outer
radius of $B^-$ is less than the inner radius of $B^+$.
Suppose we could
find a set of coefficients such that
$f_\beta({\bf r})$ satisfied
\begin{equation} \label{firstconstraint}
f_\beta({\bf r}) \le 0, \;\; {\bf r} \in B^-,
\end{equation}
\begin{equation}
f_\beta({\bf r}) \ge 0, \;\; {\bf r} \in B^+,
\end{equation}
\begin{equation}
f_\beta({\bf r}) = 0, \;\; {\bf r} \in B^0,
\end{equation}
\begin{equation}
\int_{B^-} f_\beta({\bf r}) d {\bf r} = -1 ,
\end{equation}
and
\begin{equation} \label{lastconstraint}
\int_{B^+} f_\beta({\bf r}) d {\bf r} = 1 .
\end{equation}
Then $\int_B f({\bf r}) \Omega({\bf r}) d {\bf r} < 0$ implies that in some subset of
$B^-$, $\Omega$ is larger than it is anywhere on $B^+$, and hence that the
rotation rate must somewhere increase with depth. 
The linear combination of
the data
$\sum_{j=1}^n
\delta_j$ is an unbiased estimate of $\int_B f({\bf r}) \Omega({\bf r}) d {\bf
r}$ with
variance
$\sum_j
\beta_j^2 \sigma_j^2$, and we can readily construct a confidence interval for
$\int_B f({\bf r}) \Omega({\bf r}) d {\bf r}$.  
It remains to find the coefficients $\{\beta_j\}$.
If more than one set of coefficients satisfies the restrictions (\ref{firstconstraint} - \ref{lastconstraint}), we would
like the one that minimizes $\sum_j
\beta_j^2 \sigma_j^2$, as the resulting confidence interval would be shortest.

We are thus led to the optimization problem:
\begin{equation}
\min \sum_j \beta_j^2 \sigma_j^2 
\end{equation}
subject to (\ref{firstconstraint} - \ref{lastconstraint}), which is a
semi-infinite quadratic program.
This problem can be solved using the algorithm given by Stark and Parker
(Comp. Stat., 1995) using a strategy of successive refinement of a mesh 
of constraints.  The linear inequality constraints can be imposed at a grid 
of points, and the corresponding finite-dimensional quadratic programming 
problem solved.
One may search the solution to see whether it violates the constraints at
points not in the grid.  If so, new points can be added to the grid near
the violation, and the resulting finite-dimensional problem solved.
This can be repeated until the constraints are satisfied off the grid as
well as on.  Solutions at the successive levels should be similar, so the
cost of additional iterations is expected to be modest. 
The dimension of the problem can be further reduced by using only a subset of
the data mappings, and adding new ones using a greedy criterion that checks how
nearly collinear the violation of the constraint is with the function to be
added. 
 
Sekii, Genovese, Gough and Stark (ESA, 1995) and Gough, Sekii, and Stark
(Ap.J., 1995, submitted) show using averages of the Big Bear Solar Observatory 
(BBSO) data for the years 1986 and 1988-1990 that it is possible to find 
linear combinations with isolated sign changes in colatitude, using the 
SOLA method described by Pijpers and Thompson (Astron. Astroph., 1992, 1994)
They used a fairly crude construction of the kernels that I propose to improve 
using the optimization formulation presented above.  Applying this approach to 
larger data sets from SOHO-SOI should yield precise tests of hypotheses about 
the Sun's differential rotation.

The optimization problem is very memory-intensive, as the sign constraints need
to be imposed at a relatively dense grid in the solar interior.
The problem that must be solved is a sequence of quadratic programs, with a
large number of variables (up to 10,000) and many more constraints.
Many of the constraints will be redundant, but  {\em a priori} it is hard to
tell which.  Finding an efficient means of solving this problem
will be an important part of the work.
At present, one promising approach is to send subsets of the data mappings to
many different processors on the Statistical Computing Facility (SCF) network,
and to combine and winnow subsets  using a measure of merit that includes the
degree to which the constraints are violated and the variance of the resulting
functional
$f_\beta$. As each processor completes its subproblem, it would be allocated a 
new subset
of the mappings constructed from the most successful so far.
This strategy should lend itself well to distributed computing, as each
processor has a lot of work to do between times it communicates with others;
local experience suggests that inter-process communication is the most severe
bottleneck in obtaining linear scaling of speed with the number of machines
employed.
After a screening process like this one reduces the number of potentially
useful data mappings to a reasonable number, those that remain can be used in a
larger optimization problem on a single machine with a large amount of memory.
In this way I hope to find a nearly optimal solution.

3) Physical Constraints in Rotation Estimates.

There are a number of physical constraints that can be used to improve the
reliability of rotation estimates.  For example, the rotational kinetic energy
is a quadratic functional of the rotation, and clearly must be less than the
rest mass of the Sun.  While such a constraint might seem hopelessly loose,
Backus (Geophys. J., 1989) and Stark (Geophys. J. Int., 1992) show that a 
similarly prima facie useless constraint, that the energy stored in Earth's 
magnetic
field is less than Earth's rest mass, is useful in estimating spherical harmonic
components of Earth's magnetic field.
Other physical constraints one might use include the observed oblateness
(another quadratic functional), the constraint that rotation does not cause
particle velocities to exceed the speed of light, or the constraint that
matter is gravitationally bound.  The last two constraints are pointwise
linear inequalities. Stark (J. Geophys. Res., 1992) shows how to use both 
quadratic and linear inequalities to reduce the uncertainties in inverse 
problems.  With Christopher Genovese (Carnegie-Mellon University) and Douglas 
Gough (University of Cambridge) I have begun to explore these constraints to 
test hypotheses about the distribution of angular velocity within the Sun by 
seeking models that violate various physical hypotheses while
satisfying the constraint and fitting the data.  We have developed theory that
allows us to take into account approximation errors from representing the solar
rotation as a finite linear combination of known functions; the theory appears
to extend to a wide variety of hypothesis-testing problems in
infinite-dimensional spaces.  I propose to continue this work and perform
the computations using the oblateness constraints to make inferences from 
SOHO-SOI data. 

Applying these constraints will requires special efficient
algorithms for nonsmooth convex optimization
problems in many variables, algorithms that exploit the special
structure of the problems.
A promising approach is to use a multistage algorithm that starts with
subsets of the data, keeping the dimension of the problem low,
then successively refines the solution until all data are
included.
Tests of different hypotheses could be distributed efficiently over the SCF
network, as the resulting problems share the same data but do not need
inter-process communication.

4) Asymptotic Inversion for Radial Soundspeed Structure

Let $\omega_{kl}$ be the eigenfrequency
with radial number $k$ and angular number $l$.
It can be shown (see Brodsky and Levshin, Geophys. J. Roy. astr. Soc., 1979) 
that 
\begin{equation} \label{duvalleqn}
	\frac{\pi (k+\alpha )}{\omega_{kl}} =
		\int_{r_t}^{R} \left ( \frac{r^2}{c^2(r)} -
		\left (\frac{l+1/2}{\omega_{kl}}\right )^2 
		\right )^{1/2} r^{-1} dr + O(\omega_{kl}^{-2}),
\end{equation}
where $R$ is the solar radius, $r$ is distance from the center of Sun, 
$c(r)$ is the soundspeed at radius $r$, and $\alpha$ is the
phase shift factor that captures the effect of the
boundary condition at $r=R$.
Equation (\ref{duvalleqn}) is sometimes called the Duvall equation
(after Tom Duvall). 
To first order in $\omega_{kl}^{-1}$, this relation is 
equivalent to Abel's
equation. 
The asymptotic inversion of free-oscillation data via
Abel's equation was first done for Earth by Brodsky and
Levshin (the first paper in English is Geophys. J. Roy. astr. Soc., 1979).
The technique has also been used for the inversion of solar free-oscillation 
data by Christensen-Dalsgaard and others (Nature, 1985), and by
Brodsky and Vorontsov (IAU, 1988).

These studies invoke an analytic inversion formula that, if the data
were perfect and complete (if a continuum of normal modes were observed,
without observational errors), would find
the unique one-dimensional soundspeed profile consistent with the
observations, within the asymptotic theory. 
Unfortunately, there is only a finite
(though quite large) number of data, and the data contain observational errors.
In this circumstance an infinite number of soundspeed
profiles are consistent with the observations; which of these is found by a
particular study depends critically on the interpolation scheme used to construct a
continuous curve from the noisy observations---different interpolations can
produce radically different models of soundspeed.  The
previous studies construct individual models consistent with the observations,
but as there are infinitely many such models, properties  of any particular
model may not be shared by the actual solar  soundspeed profile.  There is no
reasonable error analysis of this inversion procedure.

Garmany (Geophys. Res. Lett., 1979) proposed a method for finding an
envelope containing all solutions to Abel's problem that fit
the data adequately.
His work was extended to spherical geometry and to
incorporate additional constraints on monotonicity of seismic
velocity with depth by Stark {\em et al.} (J. Geophys. Res., 1986), and to
different measures of misfit by Stark and Parker (J. Geophys. Res., 1987).
The theoretical underpinnings of this general approach
to inference are set out by Stark (J. Geophys. Res., 1992),
who also
addresses the problem of solving the infinite-dimensional
problem rigorously, and some extensions relevant to helioseismology.
See also Lang (IEEE-IT, 1985) and Stark (Geophys. J. Rpoy. astr. Soc., 1987).

This ``strict bounds'' approach 
uses the data, the statistics of the data errors,
the mathematical relation between the data and the unknown soundspeed 
profile, and
physical constraints such as monotonicity of the soundspeed profile, 
to define a confidence region for the true soundspeed profile 
within the 
infinite-dimensional space of all soundspeed profiles.
Properties of that set are, with a specified probability,
properties of the true soundspeed.
By solving optimization problems to find the smallest and largest 
values
``interesting'' functionals attain within the confidence region, 
we get confidence intervals for functionals of the true 
soundspeed.
In particular, one can find confidence intervals for the depth at which a certain
soundspeed occurs.
This can be done, with simultaneous coverage probability, 
at an arbitrary number of soundspeeds.
The result is that one may construct a ``confidence envelope'' that, 
with specified probability, contains the true soundspeed profile.
It is possible to do this quite rigorously, incorporating the true infinite-dimensional
uncertainty, the data errors, modeling errors such as neglected terms 
in the
asymptotics, {\em etc.}
For example, in equation (\ref{duvalleqn}), both the limit of integration and 
integrand depend on $\omega_{kl}$, which is not known 
precisely---only estimates with errors are available.
Thus we do not even know precisely what functional of the Sun 
was measured.
I call this the ``fuzzy functional'' problem; Stark (J. Geophys. Res., 1992)
shows how to incorporate this additional uncertainty to get a conservative
confidence  envelope for soundspeed.
Expressions for the higher-order terms in the asymptotic relations, including
perturbations to gravity, have been derived recently by Brodsky and
Vorontsov (Ap. J., 1993); by evaluating these expressions for
plausible solar structures we can obtain reasonable bounds on the neglected
higher-order terms. The ``strict bounds'' approach can incorporate this
additional source of uncertainty as well.

The usual misfit criterion used to define a confidence region for the
noise-free observations is based on an $l_p$ norm of the errors, with $p$ equal
to 1, 2, or $\infty$. (The chi-square, or 2-norm measure is the most
common.) This procedure is not statistically consistent: as the number of data
grows, the confidence intervals for a given functional do not shrink to zero
(Stark, Phys. Earth. Planet. Int., 1995, in prep.) D.L. Donoho (personal 
communication, 1994) has proposed
using  a wavelet-based measure of misfit instead of an
$l_p$ norm (Stark, Herron, and Matteson, Appl. Spect., 1993, use such a 
measure), and claims that the wavelet misfit measure produces confidence 
intervals whose widths shrink to zero at essentially the asymptotically 
optimal rate as the number of data grows, and that automatically adapts to 
the inherent smoothness of the object.  Work by Stark also suggests that 
simply averaging sets of observations that sample the Sun in similar ways 
can also lead to optimal rates of convergence of the confidence intervals. 
Benjamini and Stark (J. Amer. Stat. Assoc., 1996, in press) show that unusual 
geometries for confidence sets can be employed to improve tests of multiple 
hypotheses, for example, improving the ability to distinguish several effects 
from zero simultaneously.

I propose to
explore and apply these approaches to estimate radially-averaged soundspeed
within the Sun, incorporating the effects of data errors, asymptotic
approximations, and ``fuzzy functionals'' into the analysis of the uncertainty
of the results.  The new misfit measures give rise to semi-infinite linear or
quadratic programs.   Stark (J. Geophys. Res., 1992)
discusses some strategies for solving such problems computationally, and 
Stark and Parker (Comp. Stat., 1995) present an algorithm for solving the 
subproblems that arise.  The dual approach is particularly
well-suited for distributed computations on the SCF network.

An important fringe benefit of this approach, which should extend to other
inverse problems as well, is that the new misfit functionals afford a
substantial ``compression'' of the observations, while enhancing our ability to
recover solar structure.
For example, Stark (Phys. Earth Planet. Inter., 1995, in prep.) proves that 
in a particular nonparametric regression problem, it is preferable to compute 
with $n^{1/3}$ averages of the data than with the $n$ original data:  the 
derived confidence intervals are shorter, and the computational savings 
are enormous.  Hengartner and Stark (Ann. Stat., 1995) proved that a 
similar phenomenon holds when estimating probability density functions 
(using a collection of weighted averages that number only a fractional 
power of the observations gives optimally narrow confidence intervals); and 
Stark, Herron and Matteson (Appl. Spect., 1993) found that a reduction of 
3601 data to 90 weighted averages (specially chosen wavelet coefficients) 
improved the recovery of rock mineralogy from Fourier-Transform infra-red 
spectroscopic data.

Exploring such ``teleological'' data compression, where the goal is to recover
solar structure, rather than to reconstruct the original observations, is
perhaps the most important part of the planned research in the long-term.
As the number of high-quality observations of solar free-oscillations grows
through SOHO-SOI, the computational savings will become crucial if the
helioseismic community is to extract as much information as possible about the
solar interior.

5) Estimating the Surface Phase Shift $\alpha$.

In order to apply the analytic inversion or strict bounds, it is necessary to
estimate the ``surface phase shift'' $\alpha$, which is a function of frequency
$\omega$. 
The phase shift is related to boundary conditions
at the solar ``surface''---see equation \ref{duvalleqn}.
The inability to remove completely the effects of $\alpha$ from
the normal mode data limits the residual noise
level in the data used in asymptotic inversions; the ``strict bounds''
formulation can incorporate this source of uncertainty.
The ratio
$\frac{\pi(k+\alpha)}{\omega_{kl}}$ is known to be monotonic {\em a priori.}
This monotonicity can be exploited to reduce the uncertainty in estimates of
$\alpha$.
Monotonicity constraints are a special case of pointwise linear inequality
constraints, as they constrain the generalized derivative to be of one sign
everywhere; thus, studying this problem should help us in more complicated
situations such as using the constraint that matter is gravitationally bound
in estimating the helioseismic rotation.

The current method for determining $\alpha$ is to fit a
simple function to the data in such a way as to minimize the scatter of the
reduced data. The method does not provide uncertainty estimates.
The data are insensitive to shifts of the entire $\alpha$ curve, so without
additional constraints, only changes in $\alpha$ as a function of frequency are
identifiable.

Many authors have studied the monotone regression problem (see
Robertson {\em et al.} (John Wiley and Sons, 1988) for an overview; see
Mukarjee and Stern (J. Amer. Stat. Assoc., 1994) for a recent development and
references), which can be summarized as follows: We observe
$n$ data
$\delta_j$,
$j=1,
\ldots, n$, that are samples of a monotone function $f$, corrupted by additive
Gaussian noise:
\begin{equation}
\delta_j = f(t_j) + \epsilon_j, \;\;j=1,\ldots, n ,
\end{equation}
where $\{\epsilon_j \}$ are independent, identically distributed, zero-mean
Gaussian random variables with unit variance.
Various techniques have been proposed for estimating $f$ from these
observations, but none gives conservative confidence intervals for
$f$ nor a simultaneous confidence envelope for $f$.

Hengartner and Stark (Ann. Stat., 1995) develop a technique for estimating
probability density functions nonparametrically, subject to a shape
restriction on the density, for example, the constraint that the
density is monotone or has at most
$k$ modes.  Their method gives a conservative confidence envelope for the 
density
whose width converges to zero at the optimal rate at all points where the
density is continuous. The basic idea behind their method is to set up a
nonparametric confidence region for the unknown density based on a
distribution-free test, then solve optimization problems subject to the shape
restriction on the density to find the largest and smallest values any density
can have and adequately fit the data and satisfy the shape restriction.
This is, in essence, the same as the ``strict bounds'' approach outlined
above.
The optimization problems can be reduced exactly to finite-dimensional linear
programs. 

Their general approach can be adapted to the monotone regression
problem by considering a confidence set for the entirety of $f$ based on a
measure of misfit to the data 
$\delta = (\delta_1, \ldots, \delta_n)$, intersecting that set with the set of
monotone functions, and solving a constrained optimization problem to find the
largest and smallest values at the point $t_0$ among functions in that
intersection. The result is a conservative confidence interval for $f(t_0)$, and
as
$t_0$  varies, the
confidence intervals have simultaneous coverage probability.
The intervals can be interpolated conservatively using the
monotonicity constraint to obtain a confidence envelope for the entire function
$f$; {\em i.e.}, a pair of functions $\ell(t)$, $u(t)$ such that
\begin{equation}
P( [\ell(t), u(t) ] \ni f(t) \;\; \forall t ) \ge 1-\alpha.
\end{equation}
As noted above, basing the
confidence set for
$f$ on a chi-square measure of misfit to the entire data vector $\delta$
is statistically inconsistent, in the sense that the
confidence intervals do not shrink as the number of data increases.
However, by
averaging the data in adjacent bins and basing the confidence set on a norm of
the misfit to those averages, one can construct consistent simultaneous
confidence intervals. Furthermore, the optimization problems
one must solve can be reduced without approximation to a sequence of
finite-dimensional linear or quadratic programs (for the one, two, or
infinity-norm),  which can be solved efficiently and accurately in practice
(Stark and Parker, Comp. Stat., 1995).

The dependence of the number of averages on the number of data determines a
tradeoff between the variance of the averages, and the bias in the averages as
estimates of the value of $f$ at the endpoints of the intervals from the
variation of $f$ on the intervals.
This is a particular instance of ``teleological data compression'' in which an
optimal way to compress the data to sharpen inferences about $\alpha$ can be
determined explicitly and fairly easily.  I propose to apply this procedure to
estimate the surface phase shift
$\alpha$ from GONG and SOHO-SOI data as they become available.  I will use the
resulting confidence intervals for
$\alpha$ as a function of $\omega$ in the asymptotic ``strict bounds''
inversion for soundspeed as a function of radius in the Sun, incorporating the
uncertainty in the estimates of $\alpha$ into the uncertainty of the estimated
soundspeed.



Philip B. Stark\\                      
Dept. of Statistics\\              
University of California\\              
Berkeley CA 94720-3860\\
stark@stat.berkeley.edu\\
 510-642-1430\\
510-642-7892 (fax)
\end{document}