# PowerGraph: Using neural networks and principal components to determine
multivariate statistical power trade-offs

Ajinkya K Mulay, Sean Lane, Erin Hennes
Department of Psychological Sciences
Purdue University
610 Purdue Mall
West Lafayette, Indiana, 47907
{mulay, lane84<EMAIL_ADDRESS>

###### Abstract

Statistical power estimation for studies with multiple model parameters is
inherently a multivariate problem. Power for individual parameters of interest
cannot be reliably estimated univariately since correlation and variance
explained relative to one parameter will impact the power for another
parameter, all usual univariate considerations being equal. Explicit solutions
in such cases, especially for models with many parameters, are either
impractical or impossible to solve, leaving researchers to the prevailing
method of simulating power. However, the point estimates for a vector of model
parameters are uncertain, and the impact of inaccuracy is unknown. In such
cases, sensitivity analysis is recommended such that multiple combinations of
possible observable parameter vectors are simulated to understand power trade-
offs. A limitation to this approach is that it is computationally expensive to
generate sufficient sensitivity combinations to accurately map the power
trade-off function in increasingly high-dimensional spaces for the models that
social scientists estimate. This paper explores the efficient estimation and
graphing of statistical power for a study over varying model parameter
combinations. We propose a simple and generalizable machine learning inspired
solution to cut the computational cost to less than 10% of the brute force
method while providing F1 scores above 90%. We further motivate the impact of
transfer learning in learning power manifolds across varying distributions.

##  1 Introduction

Statistical power is quantitatively equal to the probability of rejecting a
null hypothesis. Thus, it plays a crucial role in finding an effect of
hypothesized interest in any study. Historically, the simplest and most direct
way to improve power is to increase the study’s sample size, as
mathematically, it has the most significant impact concerning what the
researcher can control Cohen (1992). However, there are practical
considerations to make when increasing the sample size. In the case of a rare
disease study, the total population itself might be small, difficult to
locate, and even more challenging to engage, as in early- and late-stage
neuro-genetic syndromes Button et al. (2013); Szucs & Ioannidis (2017). On the
other hand, there could be access or funding issues with increasing the number
of participants a researcher wishes to obtain.

In this article, we empirically show that tuning the model parameters (_i.e._
, weights) can dramatically change the study’s power. Thus, the change in
model parameters can push the survey from a mid-powered study to well-powered
research even for a constant sample size. For the rest of the article, we
assume that well-powered research is one where the power is higher than 0.8.
Fig. 1 demonstrates that even for a fixed sample size, the power can
considerab1y vary due to a change in model parameters. In Fig. 1(a) we compute
power gradients by considering the directional difference between clusters
derived by KMeans, where the power and the squared L2-norm of the model
parameters represent the point coordinates.

Thus, generating such a manifold, as shown in Fig. 1(a) can exceptionally aid
researchers in identifying areas of high power. Enabling access to these plots
is thus a simple way of reducing the ideal sample size without sacrificing
power. However, as described in Algorithm 1, computing power even once
requires 200-1000 simulations. Further, the manifold parameter space,
including the model weights, sample size, and additional hyperparameter
choices, is enormous. For instance, for a seven predictor model, with 11
choices per predictor and 11 choices for the $N$, we have a parameter space of
$11^{7+1}$. Computation for this vast space might even take days.

Thus, in this work, we look at cheaper alternatives to simulating the entire
power manifold. We present three different contributions

  * •

We develop a neural network-based approach to predict power over a high-
dimensional manifold for classification and regression purposes. Such an
approach provides significant performance with even 10% of the training data.

  * •

We provide four baselines for alternatively predicting power and test all our
approaches on three models-= regression, logistic and Repeated Measures ANOVA.

  * •

To leverage the transfer of knowledge from one model’s power manifold to
another model’s manifold, we first demonstrate the intuition behind such a
transfer. Next, we show that such a transfer can enable the same or better
results in every case. We organize the rest of the article as follows. We
first detail the relevant literature and highlight the background information
required to build our algorithm. We provide pseudocodes for the core and
auxiliary algorithms in this section. Next, we present our empirical results,
highlight the significant wins, and report the limitations. Finally, we also
present potential future work.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1: Impact of Model Parameters on Power: (a) Change in power over
varying model coefficients ($\beta),N=105$. We observe significant power
gradient over clusters derived by KMeans for a 3-predictor regression model
with a partial f-test. (b, c, d, e, f) For the same partial f-test and over
popular regression datasets KumarRajarshi (2015); Eso (2019); Choi (2018); Too
(2019); Par (2022) ($N=50$), we demonstrate that if we perturb the estimated
model parameters ($\hat{\beta}$) by even small amounts, the under-powered
studies can become well-powered (_i.e.,_ power $>0.8$).

##  2 Related Work

Previous articles (Bakker et al. (2012), Bakker et al. (2016), Maxwell (2004),
Cohen (1992)) have pointed out that studies are frequently underpowered and
thus lead to statistically insignificant results. Underpowered studies are
primarily a result of a lack of a formal power analysis. Furthermore, in Card
et al. (2020) the authors demonstrated that numerous NLP models are typically
underpowered. They mention that datasets with sentences less than 2000 lead to
a power of only 75% and provide best practices to perform power analyses.
Next, the work in Koehn (2004) provides ways to measure the statistical
significance of the performance improvement due to a system change in the
domain of machine translation. The authors argue that just a change in the
system metric (ex., accuracy or F1-score) is not enough, and we need to
demonstrate that the change in the metric is statistically significant. This
work further broadens the need for thorough power analysis. Finally, Koehn
(2004) takes into consideration the effects of model complexity in extracting
the maximum information from a given dataset. Under a given class of models,
the authors provide an algorithm to identify the right-fit model. The right-
fit model avoids both underfitting and overfitting and has predictive power
asymptotically close to the best model in the given model class. Such a method
can replace cross-validation and signifies the importance of improving
predictive power.

Bakker et al. (2012) suggested power simulations for studies using
Questionable Research Practices (QRPs) John et al. (2012) and lower sample
sizes with more trials. QRPs could include running multiple tests with a much
smaller sample size (underpowered), rerunning analyses after adding more
subjects, or (even selectively) removing outliers. The results show that such
practices can significantly inflate statistical power and provide misleading
evidence about the effect size while hampering the study’s reproducibility.
Bakker et al. (2012) suggests that to avoid such underpowered studies, we
should use sample sizes derived after a formal power analysis.

Baker et al. (2021) solved the power issue by providing an online tool with
power contours to demonstrate the effect of sample size (N) and trials per
participant on statistical power (k). With the results provided in the
article, it is clear that changes to N or even k can dramatically change the
power region and convert an underpowered study to an appropriately powered
one. Rast & Hofer (2014) also demonstrates the strong (inversely correlated)
impact of sample size on both the effect size and study design.

Lane & Hennes (2018) and Lane & Hennes (2019) provide a clear guide to
conducting formal power analysis based on simulation methods (rather than a
formula-based approach). Even though the simulation approach is universal, it
often requires many computational resources. The resource usage increases
exponentially with a linear increase in the number of predictors or other
factors impacting power. Thus, in this work, we highlight machine learning-
based approaches that can seriously reduce resource usage with accuracies
greater than 95%.

##  3 Proposed Work

###  3.1 Computing Power

We present a power computation algorithm in Algorithm 1 for t-tests in a
linear regression model. We can extend this algorithm to the f-test by
replacing the t-test in Line 9 with an f-test or a partial f-test. For
extensions to other models, we only need to modify the data generation process
on Line 5. Computing power for multiple model weights requires us to call the
COMPUTE-POWER function each time. Thus, the cost of computing power boils down
to the number of calls made to COMPUTE-POWER. We wish to reduce the number of
calls while still predicting power for the entire parameter space in our work.

###  3.2 Feature Engineering with PCA

Complex datasets include several features, and it often becomes necessary to
reduce data dimensionality to conserve resources or speed up training.
Furthermore, we wish to remove redundancy in features. Removing highly
correlated features can improve training speeds or data processing speeds,
reduce bias, and improve the interpretability of our dataset. A common way to
achieve these goals is to use Principal Component Analysis (PCA) Wold et al.
(1987). For brevity, we add further details of PCA to the appendix. We also
provide a simple PCA algorithm in Algorithm 6 in the appendix for
completeness.

We already know that increasing the sample size $N$ can increase the power of
a study. To visualize the impact of each parameter in the model, we take a
look at a partial f-test. The data follows the distribution of the first three
variables from Table 5. We compute the power of the partial f-test to test for
the significance of the first and third predictors.

To visualize the importance of each parameter, we draw a correlation plot in
Figure 13. We can easily verify that as expected $\beta_{1}$ and $\beta_{3}$
have higher correlation with the power (as the hypothesis depends on them).
$N$ has a high correlation as we would expect. Further, we define a new
feature scaled weight computed as $N\sigma$, wherein
$\sigma=\sqrt{\sum_{i={1}}^{k}\beta_{i}^{2}}$ ($k$ denotes the number of
predictors). We observe that the $N\sigma$ feature has the highest correlation
amongst non-PCA features. After including the principal components the
correlation substantially increases to $\geq 0.90$ for $\text{PC}_{1}$. We
believe that transforming the data with PCA can further increase the variance
while ignoring the redundant features in our dataset. We, provide exact
details about our dataset in 4.2.

###  3.3 Primary Approach: The Neural Network Approach

We employ a simple neural network as described in Table 4 trained on the new
dataset. With the flexibility of a neural network, we can train for both a
classification and a regression task. We denote this approach by POWER-NETWORK
or POWER-Neural Network (PNN)). PNN also allows for an easy extension to a
multi-class classifier. However, we do not explore this domain in this
article.

####  3.3.1 Learning Faster with Transfer Learning

The primary goal of our article is to reduce the number of calls to the power
function and thus reduce the simulation overhead for computing power. Thus, we
now describe a commonly used technique fine-tuning Houlsby et al. (2019) in
Natural Language Processing (NLP) domains. Rather than cold-starting model
training with random weights, we take the help of a previously trained model
to initialize the model weights. Thus, most layers are already trained, and
our new smaller dataset can be used to tune the neural net for our new task.
Note that when the number of input dimensions is different we remove the
unused dimensions (from the pre-trained model) by specifying zero input
columns for those features.

Formally, the transfer learning algorithm P-Transfer only deviates from PNN in
the fact that line 6 of the algorithm 4 has an additional step of initializing
the model (before training) with the weights of a larger pre-trained model.
Note that the original pre-trained model needs more features than the new
model. Finally, the extra features in the larger model are turned off by
passing zero vectors to these features. Thus, even a larger model can be used
for our smaller models.

####  3.3.2 Intuition behind Transfer Learning

The intuition behind transfer learning may be sought from this article Chinn
(2000). We first note that each model directly interacts with the linear
model. However, due to the difference in their composition, the effect size
might differ. The authors inChinn (2000) empirically demonstrate that given
the standard deviation, we can find the correlation between, say, the effect
size of the logistic regression model and the regression model. Thus, we would
have the same standard deviation and thus correlated/comparable effect sizes
for standardized predictors. Thus, we might be able to use transfer learning
to boost a model with similar traits as long as the predictors are
standardized and they both have the same hypothesis.

###  3.4 Baselines: Power Cluster, Power Label Propagation and KNeighbors
Classifier

POWERCLUSTER also referred to as P-CLUSTER) simply runs K-Means on $PC_{1}$
and $N\sigma$ for identifying two clusters with two vectors from our dataset.
We provide its pseudocode in Algorithm 2. K-Means identifies clusters of
points that minimize the intra-cluster distances. Since both of our feature
vectors are independent of power, we can cluster our data points without
computing the true power.

The intuition behind deploying P-CLUSTER is that it can quickly identify
appropriately varying power domains, as seen from Fig. 2. If we compare the
clustering to the standard definition of a high-powered study (_i.e.,_ power
is more significant than 0.8), we can segregate power with high performance.
We demonstrate the results in Figures 4(a), 4(b). We can see that P-CLUSTER
performs well arbitrarily. Thus, we can conclude that not including label
information will lead to poor performance for models.

We provide two other strong baselines for power prediction- the label
propagation algorithm from Zhu & Ghahramani (2002), and the KNeighbors
Classifier Nelli (2018). We use the Label Propagation version implemented in
Scikit-Learn Pedregosa et al. (2011). Label propagation works by first
creating a fully-connected graph between all data points wherein each edge
weight is directly dependent on the euclidean distance between the two points.
These edge weights dictate the probability of propagating a label onto another
label. We direct the reader to section 2.2 of Zhu & Ghahramani (2002) for
further details on the Label Propagation algorithm. We also refer to it as PL-
PROP.

KNeighbors classifier Nelli (2018) also provides good performance, and we thus
include it for completeness. Overall, Label Propagation performs consistently
better than the KNeighbors classifier.

Algorithm 1 Computing Power of a t-test for a Linear Regression Model

1:Input: Distributions of columns in the
dataset-$\mathcal{D}=\\{D_{X_{1}},D_{X_{2}},...,D_{X_{p}}\\}$, Sample
Size-$N$, Model Weight-$\beta\in\mathbb{R}^{p}$, Number of predictors-$p$,
sensitivity-$\alpha$ (default: 0.05), number of simulations-sims (default:
1000), Error Distribution-$\mathcal{E}$

2:Output: Power of t-test

3:procedure COMPUTE-POWER($X,k$)

4: significance $\leftarrow$ 0

5: for 1:sims do

6: $X\leftarrow$ generate $N$ data samples from

7: distribution($\mathcal{D}$)

8: $e\leftarrow$ generate error from distribution($\mathcal{E}$)

9: $y\leftarrow(X\times\beta)+e$

10: $\mathcal{M}\leftarrow$ Fit $(X,y)$ to a linear regression model

11: if p-value of t-test for model $\mathcal{M}\leq\alpha$ then

12: significance += 1

13: power $\leftarrow\frac{\textit{significance}}{\textit{sims}}$

14: return power

Algorithm 2 Unsupervised clustering of the power surface with PCA features

1:Input: Length of training set-$L$, model parameter space-$\mathcal{S}$ of
length $L$ (each parameter consists of the weight vector $\beta$ and sample
size $N$), power sensitivity-$\alpha$ (default 0.05), number of simulations-
sims (default 200), PCA variance (in %) to be retained-var, Number of
clusters-$k$

2:Output: Power Surface Graph

3:procedure POWERCLUSTER($\mathcal{S},\alpha$, sims)

4: $\mathcal{S}_{\text{PCA}}$ $\leftarrow$ PCA-FIT-TRANSFORM(S, variance=var)

5: $C_{k}$ $\leftarrow$ KMeans($\mathcal{S}_{\text{PCA}}$, num_clusters $=k$)

6: return clusters $C_{k}$ derived from KMeans

Algorithm 3 Training Data Collection

1:Input: Length of training set-$L$, model parameter space-$\mathcal{S}$ of
length $L$ (each parameter consists of the weight vector $\beta$ and sample
size $N$), power sensitivity-$\alpha$ (default 0.05), number of simulations-
sims (default 200)

2:Output: Training Set $\mathcal{S}$ with parameters and corresponding powers

3:procedure GENERATE-DATA($\mathcal{S},\alpha$, sims)

4: for parameter in $\mathcal{S}$ do

5: ($\beta,N)\leftarrow$ parameter

6: power $\leftarrow$ COMPUTE-POWER($\beta,N,\alpha,$sims)

7: $\mathcal{S}$[powers] $\leftarrow$ power

Algorithm 4 Classifying/Predicting power surface with PCA features

1:Input: Length of training set-$L$, model parameter space-$\mathcal{S}$ of
length $L$ (each parameter consists of the weight vector $\beta$ and sample
size $N$), power sensitivity-$\alpha$ (default 0.05), number of simulations-
sims (default 200), PCA variance (in %) to be retained-var

2:Output: Power Surface Graph

3:procedure POWERNETWORK($\mathcal{S},\alpha$, sims)

4: $\mathcal{S}^{*}\leftarrow$ GENERATE-DATA($\mathcal{S},\alpha$, sims)

5: $\mathcal{S}_{\text{PCA}}$ $\leftarrow$ PCA-FIT-TRANSFORM($\mathcal{S}$,
variance=var)

6: Train Neural Network Model $\mathcal{M}$ on
($\mathcal{S}^{*}\cup\mathcal{S}_{\text{PCA}}$)

7: return trained model $\mathcal{M}$

##  4 Experiments

###  4.1 Code

We test the efficacy of the baseline algorithms and 4 with three kinds of
statistical models. For simulation, we use Google Colaboratory with the
standard run-time.111Our code is available at the anonymized link
https://anonymous.4open.science/r/powergraph-4BFC. Other experiment details
are deferred to the appendix in section B.1.1.

Algorithm 5 Sampling in High-Dimensions

1:Input: no. of training points-$N_{s}$, domain of the $k$ predictors and the
sample size-$\mathcal{D}\in\mathbb{R}^{k}$, no. of local training
points-$n_{s}$, local sigma for Gaussian sampling-$\sigma_{l}$

2:Output: $\mathcal{S}$ model parameter space

3:procedure P-SAMPLER($N_{s}$)

4: $n\leftarrow 0$

5: while $n<N_{s}$ do

6: Uniformly randomly pick a centroid $\beta$ from $\mathcal{D}$

7: $\mathcal{S}_{l}\leftarrow$ Sample $n_{s}$ datapoints from a Gaussian
distribution $\mathcal{N}(\mu=\beta,\sigma=\sigma_{l})$

8: $\mathcal{S}\leftarrow\mathcal{S}\cup\mathcal{S}_{l}$

9: $n\leftarrow n+1+n_{s}$ return parameter space $\mathcal{S}$

###  4.2 Dataset

For sampling parameters we use P-SAMPLER (algorithm 5) to generate the
parameter sample space and in conjunction with algorithm 3 to compute power.
Both of these algorithms together make up the power computation black box.
After collecting the dataset, we pre-process it using PCA, as explained below.
P-SAMPLER is inspired by the fungible weights described in Waller (2008).
Next, we provide exact details about the data we use for the methods
introduced in this work. For both supervised and unsupervised tasks we use the
same dataset albeit the unsupervised one does not contain the labels or true
power information. The dataset details are included in Table 1. We assume that
our original linear model is denoted as

$\displaystyle y=x_{1}\beta_{1}+x_{2}\beta_{2}+...+x_{k}\beta_{k}+e$

wherein, $e\sim\mathcal{N}(0,1)$ is the error term, $k$ is the number of
predictors, $y$ is the model predictions, $x=\\{x_{1},...,x_{k}\\}$ and
$\beta=\\{\beta_{1},...,\beta_{k}\\}$ is our linear model parameter. Our
modified dataset includes $X=\\{\beta,N,N\sigma\\}$. After including the PCA
induced features our final dataset is
$X_{PCA}=\\{\beta,N,N\sigma,\text{PC}\\}$.

####  4.2.1 Models

Three models are considered for understanding how the power manifold is
impacted by model complexity. $x$ below signifies the number of predictors in
the model.

  * •

Linear Model (REG-x): We run an f-test to understand whether a simpler model
works better than the full model. Feature distribution listed in Table 5.

  * •

Logistic Regression (LOGIT-x): Here, we run a Wald’s test to understand how
much a reduced model impacts the predictive power. This is the only non-linear
model that we consider. Feature distribution listed in Table 6.

  * •

Repeated Measures ANOVA (RMANOVA-x): We run a within-subjects RMANOVA to
measure the impact of different factors (two) on the subjects. Feature
distribution listed in Table 5.

Throughout the paper, we use the same hypothesis $\mathcal{H}_{\mathcal{O}}$
where we test whether the model parameters $\beta_{1},\beta_{3}$ are zero (in
the form of a partial F-test or a Wald’s test). However, for testing the
robustness of our transfer learning algorithm P-TRANSFER, we use the
alternative hypothesis $\mathcal{H}^{{}^{\prime}}_{\mathcal{O}}$ where we test
whether the model parameters $\beta_{1},\beta_{7},\beta_{8}$ are zero or not
for both an F-test and a Wald’s test.

Feature | Details
---|---
$\beta$ | As provided by the expected dataset’s distribution
$N,N\sigma$ | Sample Size, Scaled Weight
$\text{PC}=\\{\text{PC}_{1},...,\text{PC}_{r}\\}$ |  $r$ Principal Components as extracted from $X$ with 99% variance
Table 1: Dataset: Features included for model training

For the 3-predictor network we have 9 total features, 13 features for the
5-predictor network, 23 for the 10-predictor network and 43 for the
20-predictor network.

###  4.3 Metrics

In our proposed algorithm, we solve two kinds of problems- binary power
classification and a regression problem with actual power values.

####  4.3.1 Classification

For binary classification, we normally simply use the he accuracy. However,
since we might have unbalanced datasets, we instead report the F1 score over
the test set. While using multiple binary classifiers we could still use the
F1-score albeit in an one-vs-rest fashion. For multi-class classification, we
are first required to run the more restrictive binary classifier, compute its
performance, and then for the second binary classifier we only compute its
performance over the remaining data points. Thus, suppose we have two binary
classifiers- $C_{1}$ (power $>0.8$ and power $\leq 0.8$) and $C_{2}$ (power
$>0.6$ and power $\leq 0.6$). We first run $C_{1}$ to find all data points
with power $>0.8$ and then run $C_{2}$ to find data points with $0.8\geq$power
$>0.6$ and power $\leq 0.6$. Note that currently in the results we only
provide separate performance metrics for $C_{1}$ and $C_{2}$.

####  4.3.2 Regression

We use a Mean Squared Error loss function for the regression problem between
the ground truth power and the simulated power values. However, power is the
probability of successfully rejecting a null hypothesis, and thus we can
consider the ground truth power $P^{g}$ and the simulated power values $P^{s}$
as probability distributions. Note that the KL Divergence between $A$ and $B$
is shown in equation 1.

$\mathcal{L}_{KL}(A||B)=\sum_{x}A(x)\text{log}\Big{(}\frac{A(x)}{B(x)}\Big{)}$
(1)

Here, $x$ denotes the feature space associated with the probability
distribution (_i.e.,_ actual and simulated power), and $A,B$ represents the
ground truth and the simulated powers themselves. However, KL divergence is
non-symmetric and not a distance measure. Thus, we use the Jensen-Shannon (JS)
divergence which is an extension of the equation. 1. The JS divergence between
two probability distributions $A$ and $B$ is given by the equation 2.

$\mathcal{L}_{JS}(A||B)=\frac{1}{2}\\{\mathcal{L}_{KL}(A||B)+\mathcal{L}_{KL}(B||A)\\}$
(2)

(a)

(b)

Figure 2: PNN-CLUSTER Performance: A simple unsupervised learning algorithm
can identify high and low power domains. (left) Results for partial f-test
performed on a 3-predictor REG model. (right) Results for partial f-test
performed on a 5-predictor REG model.

###  4.4 A note about current baselines for predicting power

From literature, previous approaches for computing a power manifold used a
brute-force approach by essentially running the black box for computing power.
Thus, given unlimited time, they could achieve 100% accuracy by just literally
computing the power manifold. Specific greedy approaches attempted to find
optimal power but not the entire manifold. Since our work aims to minimize the
effort required to compute the power manifold, we do not compare our work
directly to these brute-force approaches. Instead, we elect to create our
baselines using traditional predictive machine learning.

##  5 Results

###  5.1 P-Transfer

We evaluate P-Transfer by first pre-training a REG-20 model and then testing
it over three different scenarios-(i) Same feature distribution, Same
Hypothesis, (ii) Different feature distribution, Same Hypothesis, and (iii)
Different feature distribution, Different Hypothesis. We denote these cases by
$TF_{1},TF_{2},TF_{3}$ respectively. Note that we always use a model smaller
or equal in size to the original pre-trained model for transfer learning.

Our preliminary observations indicate that transfer learning never hurts model
performance. However, for the cases where both the models have the same
hypothesis (_i.e.,_ cases (i) and (ii)), we see significant gains in the model
performance. Table 2 summarizes these results.

Type |  Model |  P-Transfer Performance |  Original Performance (Same Amt. of Train Data) |  Original Performance (8x Amt. of Train Data)
---|---|---|---|---
$TF_{1}$ |  REG3 |  $\bm{0.9902}^{\dagger}\bm{\pm 0.0045}$ |  $0.9827$ |  $0.9900$
$TF_{1}$ |  REG5 |  $\bm{0.9804}^{\dagger}\bm{\pm 0.0064}$ |  $0.9700$ |  $0.9800$
$TF_{1}$ |  REG10 |  $\bm{0.9630\pm 0.0087}$ |  $0.9600$ |  $0.9800$
$TF_{1}$ |  RMANOVA3 |  $\bm{0.9458\pm 0.0104}$ |  $0.945$ |  $0.9600$
$TF_{1}$ |  RMANOVA5 |  $\bm{0.9650\pm 0.0084}$ |  $0.9600$ |  $0.9700$
$TF_{1}$ |  RMANOVA10 |  $\bm{0.9465\pm 0.0103}$ |  $0.9200$ |  $0.9800$
$TF_{1}$ |  RMANOVA20 |  $0.8513\pm 0.0164$ |  $\bm{0.8600}$ |  $0.9300$
$TF_{2}$ |  LOGIT3 |  $\bm{0.9717}^{\dagger}\bm{\pm 0.0076}$ |  $0.9300$ |  $0.9700$
$TF_{2}$ |  LOGIT5 |  $\bm{0.9524\pm 0.0098}$ |  $0.9200$ |  $0.9600$
$TF_{2}$ |  LOGIT3 |  $\bm{0.9361\pm 0.0112}$ |  $0.8900$ |  $0.9700$
$TF_{2}$ |  LOGIT3 |  $\bm{0.9099\pm 0.0132}$ |  $0.8600$ |  $0.9300$
$TF_{2}$ |  REG10 |  ${0.9658\pm 0.0178}$ |  $\bm{0.9659}$ |  $0.9944$
$TF_{2}$ |  RMANOVA10 |  $\bm{0.9260\pm 0.0256}$ |  $0.8678$ |  $1.0000$
$TF_{3}$ |  RMANOVA20 |  $\bm{0.815\pm 0.0380}$ |  $\bm{0.815}$ |  $0.815$
Table 2: P-Transfer Test Performance: Models with the same hypothesis may
provide considerable performance boost. Any other scenario does not seem to
hurt transfer learning performance.$\dagger$-indicates that P-Transfer beat
even the 8x training data performance. The confidence ranges indicate a 95%
confidence interval.

Thus, for comparison to PNN, we use our baselines -PL-Rand, P-CLUSTER, PL-PROP
and the PK-neighbors approach. Note that to accommodate the various trends in
our model, we only report the best performing baseline (over all four
baselines) in the plots. For brevity and due to many baseline performance
combinations, we skip the point-wise results of the baselines. However, the
code provided hosts with all the information required to recreate these
baselines.

We capture three primary trends in our experiments, the performance over the
change in model complexity, the number of predictors, and the amount of
training data used for predicting the power manifold.

The trend of increased performance with increased training data usage is
reported via Figures 4, 3. Interestingly, even at lower training data
proportions, we receive good performance from the classifiers. From the exact
figures, we can see that predictions struggle over the RMANOVA model while
they are much easier for the REG model and moderate for the LOGIT model. A
similar trend can be observed for regression tasks as seen in the Figs. 7, 8,
9 in the appendix. Furthermore, with a change in the number of predictors as
seen in Figs. 10, 11 and 12, the model performance deteriorates for change in
the number of predictors for REG and LOGIT. However, for the RMANOVA model, we
see an arbitrary performance plot indicating our sub-optimal model choice
might be the reason. Here sub-optimal indicates the model derived without
using the validation set. It seems that for both regression and classification
tasks, PNN struggles the most for the RMANOVA model family. However, given
that we see this is only the models with low access to training data, we can
say that these classifiers are biased.

Finally, we can observe the number of calls v/s the performance in Fig. 15.
Clearly, for multiple models using fewer data points can still provide high
utility.

Even though we do not fully report the baselines, we would like to highlight a
few notable observations. As expected, P-RAND performs the worst, followed by
P-CLUSTER. Depending on the simulated data’s distribution P-CLUSTER
surprisingly beats both PK-neighbors and PL-PROP when fewer training data
points are used. Amongst, PL-PROP and PK-neighbors, the PL-PROP approach is
the more robust one and it is the closest one to PNN.

(a)

(b)

Figure 3: PNN Classification Performance: Comparison of model performance over
change in proportion of training data used (increasing from left to right).
(left) 3-predictor performance, and (right) 5-predictor network performance.

(a)

(b)

Figure 4: PNN Classification Performance: Comparison of model performance over
change in proportion of training data used (increasing from left to right).
(left) 10-predictor performance, and (right) 20-predictor network performance.

##  6 Limitations and Future Work

This work primarily investigates the utility of training a neural network for
predicting statistical power over a manifold. Furthermore, we also provide
baselines- Random (P-Rand), K Neighbors Classifier (PK-Neighbors), Label
Propagation (PL-Prop), and an unsupervised clustering method (P-CLUSTER).

Note that each of the baselines fails to generalize for complex models.

  * •

P-RAND: The random mechanism works as advertised and provides the least
performance throughout.

  * •

P-CLUSTER: This approach can be compelling for visualizing the power manifold
in 2 dimensions space effectively. However, we note that its performance is
inconsistent, and the linear classifier might not always be able to predict
exactly between high and low power regions. Further, the lack of label
information makes this approach useful for pre-processing but not for
predicting power classes. Alternatively, we might be able to identify and
eliminate low powered regions using this approach, and we again propose this
as an open problem.

  * •

PK-neighbors: This approach is a supervised clustering alternative to
P-CLUSTER. However, it does not necessarily perform consistently. We, however,
note that we have not fine-tuned this approach for the number of neighbors. We
leave this investigation as part of our future work.

  * •

PL-PROP: The most promising baseline is PL-PROP with its rbf kernel. PL-PROP
is better than the rest but fails to generalize for complex models or more
predictors. Furthermore, it requires much more data compared to the neural
network. We, however, have not used any guidance for tuning this algorithm’s
parameters, and we leave this aspect for future work.

  * •

Model Limitations: PNN provides excellent empirical performance and
generalizes well over increasing model complexity while providing competitive
performance even with only 10% of the data. We note that PNN is not yet
trained using methods like early stopping or an adaptive learning rate. We
would pursue these optimizations in subsequent work. Particularly we do not
use a validation dataset to improve the model generalization, and we set 500
epochs as an arbitrary stopping point. We have seen at least one instance of
such an approach leading to a sub-optimal model choice, thus leaving a gap in
the empirical performance.

  * •

Dataset Limitations: Our most complex model currently is the RMANOVA with two
within-subjects factors. To further understand the generalization error, we
hope to extend our approach to higher complexity models, including other non-
linear models. Further, we would like to explore multi-class classifications
to provide the user with more options in terms of choosing the power manifold.

  * •

Theoretical Limitation: Finally, the current implementations do not provide
any formal theoretical guarantees for the convergence of any of these methods.
We do provide intuition about why our transfer learning methods work. However,
we would like to explore how different distributions with the same hypothesis
impact the power manifold in greater detail.

##  7 Conclusion

In this work, we show that PCA-derived features are beneficial for exploring
the high-dimensional manifolds of the power surface. We provide multiple
algorithms that provide reasonable alternatives to the original brute force
method of a power analysis. We show that using a simple, fully connected
neural network, we can generalize across complex and even non-linear models to
consistently predict power. We showcase that even using 10% training data can
lead to high prediction accuracy for regression and classification tasks.
Finally, with transfer learning, we can learn a single model and boost the
performance of other models again with only 10% training data.

## References

  * Eso (2019) 2019\.  URL https://www.rdocumentation.org/packages/datasets/versions/3.6.2/topics/esoph.
  * Too (2019) 2019\.  URL https://www.rdocumentation.org/packages/datasets/versions/3.6.2/topics/ToothGrowth.
  * Par (2022) 2022\.  URL https://archive.ics.uci.edu/ml/datasets/parkinsons+telemonitoring.
  * Akiba et al. (2019) Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama.  Optuna: A next-generation hyperparameter optimization framework.  In _Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining_, pp. 2623–2631, 2019.
  * Baker et al. (2021) Daniel H Baker, Greta Vilidaite, Freya A Lygo, Anika K Smith, Tessa R Flack, André D Gouws, and Timothy J Andrews.  Power contours: Optimising sample size and precision in experimental psychology and human neuroscience.  _Psychological Methods_ , 26(3):295, 2021.
  * Bakker et al. (2012) Marjan Bakker, Annette Van Dijk, and Jelte M Wicherts.  The rules of the game called psychological science.  _Perspectives on Psychological Science_ , 7(6):543–554, 2012.
  * Bakker et al. (2016) Marjan Bakker, Chris HJ Hartgerink, Jelte M Wicherts, and Han LJ van der Maas.  Researchers’ intuitions about power in psychological research.  _Psychological science_ , 27(8):1069–1077, 2016\.
  * Button et al. (2013) Katherine S Button, John PA Ioannidis, Claire Mokrysz, Brian A Nosek, Jonathan Flint, Emma SJ Robinson, and Marcus R Munafò.  Power failure: why small sample size undermines the reliability of neuroscience.  _Nature reviews neuroscience_ , 14(5):365–376, 2013.
  * Card et al. (2020) Dallas Card, Peter Henderson, Urvashi Khandelwal, Robin Jia, Kyle Mahowald, and Dan Jurafsky.  With little power comes great responsibility.  _arXiv preprint arXiv:2010.06595_ , 2020.
  * Chinn (2000) Susan Chinn.  A simple method for converting an odds ratio to effect size for use in meta-analysis.  _Statistics in medicine_ , 19(22):3127–3131, 2000\.
  * Choi (2018) Miri Choi.  Medical cost personal datasets, 2018.  URL https://www.kaggle.com/mirichoi0218/insurance.
  * Cohen (1992) Jacob Cohen.  Things i have learned (so far).  In _Annual Convention of the American Psychological Association, 98th, Aug, 1990, Boston, MA, US; Presented at the aforementioned conference._ American Psychological Association, 1992.
  * Houlsby et al. (2019) Neil Houlsby, Andrei Giurgiu, Stanislaw Jastrzebski, Bruna Morrone, Quentin De Laroussilhe, Andrea Gesmundo, Mona Attariyan, and Sylvain Gelly.  Parameter-efficient transfer learning for nlp.  In _International Conference on Machine Learning_ , pp. 2790–2799. PMLR, 2019.
  * John et al. (2012) Leslie K John, George Loewenstein, and Drazen Prelec.  Measuring the prevalence of questionable research practices with incentives for truth telling.  _Psychological science_ , 23(5):524–532, 2012\.
  * Koehn (2004) Philipp Koehn.  Statistical significance tests for machine translation evaluation.  In _Proceedings of the 2004 conference on empirical methods in natural language processing_ , pp. 388–395, 2004.
  * KumarRajarshi (2015) KumarRajarshi.  Life expectancy (who), 2015.  URL https://www.kaggle.com/kumarajarshi/life-expectancy-who?select=Life+Expectancy+Data.csv.
  * Lane & Hennes (2018) Sean P Lane and Erin P Hennes.  Power struggles: Estimating sample size for multilevel relationships research.  _Journal of Social and Personal Relationships_ , 35(1):7–31, 2018.
  * Lane & Hennes (2019) Sean P Lane and Erin P Hennes.  Conducting sensitivity analyses to identify and buffer power vulnerabilities in studies examining substance use over time.  _Addictive behaviors_ , 94:117–123, 2019.
  * Maxwell (2004) Scott E Maxwell.  The persistence of underpowered studies in psychological research: causes, consequences, and remedies.  _Psychological methods_ , 9(2):147, 2004.
  * Nelli (2018) Fabio Nelli.  Machine learning with scikit-learn.  In _Python Data Analytics_ , pp. 313–347. Springer, 2018.
  * Pedregosa et al. (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011\.
  * Rast & Hofer (2014) Philippe Rast and Scott M Hofer.  Longitudinal design considerations to optimize power to detect variances and covariances among rates of change: simulation results based on actual longitudinal studies.  _Psychological methods_ , 19(1):133, 2014.
  * Szucs & Ioannidis (2017) Denes Szucs and John PA Ioannidis.  Empirical assessment of published effect sizes and power in the recent cognitive neuroscience and psychology literature.  _PLoS biology_ , 15(3):e2000797, 2017.
  * Waller (2008) Niels G Waller.  Fungible weights in multiple regression.  _Psychometrika_ , 73(4):691–703, 2008.
  * Wold et al. (1987) Svante Wold, Kim Esbensen, and Paul Geladi.  Principal component analysis.  _Chemometrics and intelligent laboratory systems_ , 2(1-3):37–52, 1987.
  * Zhu & Ghahramani (2002) Xiaojin Zhu and Zoubin Ghahramani.  Learning from labeled and unlabeled data with label propagation.  2002\.

##  Appendix A Appendix

##  Appendix B PNN Performance

We demonstrate additional results about the model’s performance across changes
in the number of predictors, for different classification boundaries (0.6
instead of 0.8), and the performance for regression tasks.

We first consider changes in the classification boundary from 0.8 to 0.6 to
observe models with lower power as potentially interesting. We observe from
Figs. 5 and 6.

(a)

(b)

Figure 5: PNN Classification Performance: Comparison of model performance over
change in proportion of training data used (increasing from left to right).
Results for the 0.6 classification boundary for (left) 3 predictors , and
(right) 5 predictors.

(a)

(b)

Figure 6: PNN Classification Performance: Comparison of model performance over
change in proportion of training data used (increasing from left to right).
Results for the 0.6 classification boundary for (left) 10 predictors , and
(right) 20 predictors.

For all regression tasks in Figs. 7, 8, 9 (even over the change in the
classification boundary from 0.8 to 0.6) we can observe that the JS divergence
over the test set is always better than the random baseline by at least 10x
and going up to even 200x performance.

The performance over the change in the number of predictors as seen in Figs.
10, 11 and 12 shows a downward trend for the REG and LOGIT models and so the
classifier consistently provides poor performance with an increase in model
complexity due to the number of predictors. This can be explained by the
massive increase in the dimensionality of the space. However, for RMANOVA we
do not see such a trend at all, and rather the models can have higher
performance for even the complex models. We believe this is a current
limitation due to sub-optimal model choice (not driven by the validation set).

(a)

(b)

Figure 7: PNN REG Regression Performance: Regression task is a super set of
the classification task and thus inherently much harder. We report the
performance of this task over change in the number of training points.
Interestingly, even with just 10% of the data the neural network performs
close enough to the neural network using 8x training data. Results for the REG
model family (left) (classification boundary 0.8) and (right) (classification
boundary 0.6).

(a)

(b)

Figure 8: PNN RMANOVA Regression Performance: Regression task is a super set
of the classification task and thus inherently much harder. We report the
performance of this task over change in the number of training points.
Interestingly, even with just 10% of the data the neural network performs
close enough to the neural network using 8x training data. Results for the
RMANOVA model family (left) (classification boundary 0.8) and (right)
(classification boundary 0.6).

(a)

(b)

Figure 9: PNN LOGIT Regression Performance: Regression task is a super set of
the classification task and thus inherently much harder. We report the
performance of this task over change in the number of training points.
Interestingly, even with just 10% of the data the neural network performs
close enough to the neural network using 8x training data. Results for the
LOGIT model family (left) (classification boundary 0.8) and (right)
(classification boundary 0.6).

(a)

(b)

Figure 10: PNN REG Classification Performance: We demonstrate how the change
in model complexity due to increase in predictors affects PNN’s performance.
(left) (classification boundary 0.8) and (right) (classification boundary
0.6).

(a)

(b)

Figure 11: PNN RMANOVA Classification Performance: We demonstrate how the
change in model complexity due to increase in predictors affects PNN’s
performance. (left) (classification boundary 0.8) and (right) (classification
boundary 0.6).

(a)

(b)

Figure 12: PNN LOGIT Classification Performance: We demonstrate how the change
in model complexity due to increase in predictors affects PNN’s performance.
(left) (classification boundary 0.8) and (right) (classification boundary
0.6).

###  B.1 Parameter Tuning

Tuning the algorithms used in this work is critical to improving the
performance of the Power Network and Power Cluster. It turns out that choosing
the variance captured by PCA directly affects the Power Network. Moreover, PCA
variance and the number of clusters selected in Power Clustering can improve
its performance.

####  B.1.1 Experiment Details

We have standardized the parameter distribution for data collection to provide
a consistent comparison. These parameters can be found in Table 3. For tuning
the neural network’s learning rate, we use the library Optuna (Akiba et al.
(2019)).

Experimentation Details | Values
---|---
No. of power simulations for each record | 1000
Total Datapoints per model | 2000
No. of training epochs | 500
Training Data Split | 10, 20, 40, 60, 80 (%)
Power Classification Boundary | 0.6, 0.8
Optimizer | Adam
Model Performance confidence interval | 95%
Table 3: Experimental Parameters

####  B.1.2 PCA

Usually, increasing the sample size improves power. Furthermore, our feature
engineered dimension denoted by $N\sigma$ _i.e.,_(scaled weight) also affects
the power significantly. We can observe the correlation values concerning
power in Fig. 13. Thus, we expect the variance covered in PCA transformed data
to require at least two dimensions.

We consistently observe excellent performance after selecting 99% variance
from PCA. Furthermore, we concatenate these new features to the original
matrix and use this to train the Neural Network. As we see from figure 13, the
impact of the PCA features is sizeable.

Figure 13: Impact of PCA: Correlation of parameters with true power. With
feature engineering and PCA we are able to obtain correlation $>0.9$ for a
linear regression model with 3 parameters.

####  B.1.3 Choice of clusters in Power Cluster

Depending on how the data is clustered, the Choice of several clusters can
directly impact the classification performance. However, empirically this
method is unreliable and requires a visual inspection of the dataset to be
helpful.

Figure 14: P-CLUSTER Performance (multi-clustering): P-CLUSTER can identify
multiple zones in the power domain quite efficiently. (left) displays 5
clusters and, (right) displays 10 clusters.

####  B.1.4 Other Considerations

We employ a few techniques during neural network training with increasingly
complex models and parameter spaces.

  * •

The training is sensitive to the Choice of the learning rate. Tuning needs to
be achieved after using a hyperparameter tuner like Optuna Akiba et al.
(2019).

  * •

Training tends to fluctuate a lot, and a higher number of epochs combined with
a validation set and early stopping seem to help.

####  B.1.5 Neural Network Architecture

Layer | Parameters
---|---
Dense (input) | 64 units + ReLU
Dense (hidden) | 32 units + ReLU
Dense (output) | 1 unit + Sigmoid
Table 4: PNN Architecture: We report the fully connected layers. The input
dimensions depend on the number of predictors and additional PCA features.

####  B.1.6 List of Model Parameters

$D_{X_{1}}$ | $D_{X_{2}}$ | $D_{X_{3}}$ | $D_{X_{4}}$ | $D_{X_{5}}$
---|---|---|---|---
CAT[-1, 1] | $\mathcal{N}(0,1)$ | $X_{1}\times X_{2}$ | $\mathcal{N}(0,2)$ | $X_{4}\times X_{2}$
$D_{X_{6}}$ | $D_{X_{7}}$ | $D_{X_{8}}$ | $D_{X_{9}}$ | $D_{X_{10}}$
CAT[0, 1, 2] | $\mathcal{N}(0,2)$ | $X_{6}\times X_{7}$ | $\mathcal{N}(0,1)$ | $X_{2}\times X_{6}$
$D_{X_{11}}$ | $D_{X_{12}}$ | $D_{X_{13}}$ | $D_{X_{14}}$ | $D_{X_{15}}$
$\mathcal{N}(0,3)$ | $\mathcal{N}(0,1)$ | $X_{1}\times X_{11}$ | $\mathcal{N}(0,2)$ | $X_{11}\times X_{12}$
$D_{X_{16}}$ | $D_{X_{17}}$ | $D_{X_{18}}$ | $D_{X_{19}}$ | $D_{X_{20}}$
$\mathcal{N}(0,2)$ | $\mathcal{N}(0,2)$ | $X_{11}\times X_{14}$ | $\mathcal{N}(0,1)$ | $X_{6}\times X_{16}$

Table 5: $D_{\mathcal{O}}$: Distribution of model features. For predictors
$k<20$, we only select the first $k$ predictors. Note that CAT refers to a
categorical variable while $\mathcal{N(\mu,\sigma)}$ is the Normal
Distribution with mean $\mu$ and standard deviation $\sigma$.

$D_{X_{1}}$ | $D_{X_{2}}$ | $D_{X_{3}}$ | $D_{X_{4}}$ | $D_{X_{5}}$
---|---|---|---|---
$\mathcal{N}(0,1)$ | CAT[-1, 1] | $\mathcal{N}(0,1)$ | $\mathcal{N}(0,1)$ | $X_{2}\times X_{1}$
$D_{X_{6}}$ | $D_{X_{7}}$ | $D_{X_{8}}$ | $D_{X_{9}}$ | $D_{X_{10}}$
$X_{2}\times X_{3}$ | $X_{2}\times X_{4}$ | $\mathcal{N}(0,1)$ | $\mathcal{N}(0,1)$ | $X_{2}\times X_{8}$
$D_{X_{11}}$ | $D_{X_{12}}$ | $D_{X_{13}}$ | $D_{X_{14}}$ | $D_{X_{15}}$
$\mathcal{N}(0,3)$ | $\mathcal{N}(0,1)$ | $X_{1}\times X_{11}$ | $\mathcal{N}(0,2)$ | $X_{11}\times X_{12}$
$D_{X_{16}}$ | $D_{X_{17}}$ | $D_{X_{18}}$ | $D_{X_{19}}$ | $D_{X_{20}}$
$\mathcal{N}(0,2)$ | $\mathcal{N}(0,2)$ | $X_{11}\times X_{14}$ | $\mathcal{N}(0,1)$ | $X_{6}\times X_{16}$

Table 6: $D_{\mathcal{A}}$: Alternate Distribution of model features (used for
non-linear models and testing purposes). For predictors $k<20$, we only select
the first $k$ predictors. Note that CAT refers to a categorical variable while
$\mathcal{N(\mu,\sigma)}$ is the Normal Distribution with mean $\mu$ and
standard deviation $\sigma$.

##  Appendix C Computational Costs

Figure 15: Cost of Power Computation: We report the calls required to compute
the power manifold over a subset of the parameter space. Visualizing from left
to right, we can see that we can yield high-performing classifiers even with a
fraction of the training data.

###  C.1 Details about Principal Component Analysis

PCA works by finding the direction with maximal variance. Next, it finds the
direction with the maximal variance such that this direction is uncorrelated
with the previous direction. We continue in this fashion so that any pair of
directions are uncorrelated. We then reorient our dataset with these new
components to compute our new dataset. Suppose our original dataset is
$X\in\mathbb{R}^{N\times p}$ and $v\in\mathbb{R}^{p\times k}$ are the derived
Principal Components (PCs) where $N$ is the total samples in $X$, $p=$dim($X$)
and $k$ is the number of PCs to be selected. Then the new dataset is given by,

$X^{{}^{\prime}}=X\times v.$

Algorithm 6 Transform the dataset into the first $k$ Principal Components

1:Input: $X\in\mathbb{R}^{N\times p}$ dataset, number of samples in $X$-$N$,
data dimension-$p$, number of components to be included with PCA-$k$ ($k<p$)

2:Output: Transformed PCA dataset

3:procedure PCA-FIT-TRANSFORM($X,k$)

4: for $X_{i}$ in columns($X$) do

5: $\mu_{X_{i}}$ $\leftarrow$ MEAN($X_{i}$)

6: $\sigma_{X_{i}}$ $\leftarrow$ STANDARD-DEVIATION($X_{i}$)

7: $X_{i}$ $\leftarrow$ $\frac{(X_{i}-\mu_{X_{i}})}{\sigma_{X_{i}}}$

8: $Y\leftarrow X.T\times X$

9: $V_{k}\leftarrow$ FIRST-k-EIGENVECTORS($Y$)

10: return $X\times V_{k}$


# Evidence of solar-cycle related structural changes in the solar convection
zone

Sarbani Basu Department of Astronomy, Yale University, New Haven, CT 06520,
USA<EMAIL_ADDRESS>
###### Abstract

While it has been relatively easy to determine solar-cycle related changes in
solar dynamics, determining changes in structure in the deeper layers of the
Sun has proved to be difficult. By using helioseismic data obtained over two
solar cycles, and sacrificing resolution in favour of lower uncertainties, we
show that there are significant changes in the solar convection zone, and
perhaps even below it. Using MDI data, we find a relative squared sound-speed
difference of $(2.56\pm 0.71)\times 10^{-5}$ at the convection-zone base
between the maximum of solar Cycle 23 and the minimum between Cycles 23 and
24. The squared sound-speed difference for the maximum of Cycle 24 obtained
with HMI data is $(1.95\pm 0.69)\times 10^{-5}$. GONG data support these
results. We also find that the sound speed in the solar convection zone
decreases compared to the sound speed below it as the Sun becomes more active.
We find evidence of changes in the radial derivative of the sound-speed
difference between the solar minimum and other epochs at the base of the
convection zone implying possible small changes in the position of the
convection-zone base, however, the results are too noisy to make any
definitive estimates of the change.

##  1 Introduction

Solar frequencies, and frequency splittings, change with solar activity. This
has been known for decades (Woodard & Noyes, 1985; Palle et al., 1989;
Elsworth et al., 1990; Libbrecht & Woodard, 1990) and is one of the earliest,
and most robust, findings of helioseismology. The frequency changes were shown
to be correlated with the average line-of-sight magnetic field over the solar
surface (Woodard et al., 1991); this was confirmed later using much longer
data sets (e.g. Howe et al., 2002). The correlation between frequency shifts
and activity is also seen for other solar activity proxies (Jain et al., 2009;
Howe et al., 2018a, etc.).

The solar-cycle dependent changes in solar oscillation frequencies have been
successfully used to determine activity-related changes in solar dynamics (see
review by Howe 2009 for Cycle 23 results, subsequent results can be found in
Antia & Basu 2013; Howe et al. 2013, 2018b; Kosovichev & Pipin 2019; Basu &
Antia 2019, etc.) While changes in solar dynamics have been quite easy to
detect, detecting changes in solar structure has been very difficult.

Solar dynamics, and changes thereof, are determined from frequency splittings,
while the structure is determined from the frequencies themselves. Solar-cycle
dependent frequency changes are predominantly a function of frequency, with
the changes increasing with frequency. It is reminiscent of the “surface
term,” i.e., changes or differences in frequencies caused by near-surface
changes or differences. The apparent lack of a degree-dependence in the
frequency shifts has led to the view that solar-cycle related changes in the
Sun are confined to a thin layer close to, or even above, the surface
(Libbrecht & Woodard, 1990; Goldreich et al., 1991; Evans & Roberts, 1992;
Nishizawa & Shibahashi, 1995, etc.). Other studies seemed to confirm this
(Howe et al., 1999; Basu & Antia, 2000, etc.). Despite this, there have been
attempts to look for changes deeper inside the Sun.

The solar dynamo is believed to be situated in the tachocline which is more or
less coincident with the base of the solar convection zone (Kosovichev, 1996;
Basu, 1997; Charbonneau et al., 1999; Basu & Antia, 2019). There has been an
expectation that the large magnetic fields generated there will be visible as
changes in solar structure close to the convection-zone base. Consequently,
there have been attempts, with various degrees of success in detecting such
changes (Dziembowski et al., 2000; Basu, 2002; Eff-Darwich et al., 2002; Chou
& Serebryanskiy, 2005, etc.). Baldner & Basu (2008) did a principal component
analysis of the helioseismic data available until May 2007 to find a change in
sound speed at the base of the convection zone; the results at other radii
were not statistically significant. None of the other efforts has shown a
clear, unambiguous pattern of changes that are correlated with solar activity.
Attempts to determine changes in relatively shallow layers were, however,
successful. Baldner & Basu (2008) showed that there was evidence of changes in
the outer 10–15% of the Sun, but at low latitudes, while Basu et al. (2007)
found evidence of change even at high latitudes, but their work was only
sensitive to the outer 2% of the Sun. Basu & Mandel (2004) examining
signatures caused by the acoustic glitch in the Heii ionization zones detected
a solar cycle-dependent change; this was subsequently confirmed by Watson &
Basu (2020) using solar frequencies for Cycles 23 and 24, but while the result
was obtained with data sensitive to changes at all latitudes, the change was
localized to a shallow layer.

Watson & Basu (2020) showed that the use of two solar cycles worth of data
helps in demonstrating whether or not there is a solar-cycle dependent change
when the changes are small. In this paper, we use helioseismic data for two
solar cycles to examine whether or not there are changes in the deeper layers
of the Sun. We do a differential analysis, i.e., we analyze frequency
differences between different epochs of the two solar cycles and the solar
minimum between Cycles 23 and 24. We perform solar structure inversions,
however, unlike usual inversions, we sacrifice radial resolution in order to
reduce uncertainties in the results.

The rest of the paper is organized as follows: We describe the data in Section
2, the frequency differences are examined in Section 3, the inversion method
and possible systematic errors are described in Section 4, we present our
results in Section 5. We discuss the implications of our results in Section 6.

##  2 Data used

For this work, we use solar oscillation frequencies obtained by the ground-
based Global Oscillation Network Group (GONG) (Hill et al., 1996) and the
space-based the Michelson Doppler Imager (MDI) onboard the Solar and
Heliospheric Observatory spacecraft (Scherrer et al., 1995) and the
Helioseismic and Magnetic Imager (HMI) (Scherrer et al., 2012) onboard the
Solar Dynamics Observatory.

The GONG data we use cover a period from 1995.05.05 to 2020.05.18. The data
are designated by GONG “months”, each “month” being 36 days long. Solar
oscillation frequencies and splittings of sets starting GONG Month 2 are
obtained using 108-day (i.e., 3 GONG months) time-series. There is an overlap
of 72 days between different data sets, i.e., GONG Month 2 frequencies were
obtained from data of GONG Months 1, 2 and 3, those for GONG 3 from GONG
Months 2, 3 and 4, etc. We use data for GONG Months 2 to 254. While the
overlapping sets are useful in looking at time variations, they can bias
statistics, consequently, all statistical inferences are drawn from non-
overlapping data sets. These data sets are publicly available from GONG
archives.

Data from MDI cover the period from 1996.05.01 to April 24, 2011.04.24. Solar
oscillation frequencies and splittings for these data are obtained from 72-day
time series and the sets have no overlap in time. HMI started obtaining data
on 2010.04.30 and these are also obtained from 72-day time series. We use all
MDI data, and all HMI sets up to set 10000, which has an end date of
2020.07.30. These sets are available to all researchers through the Joint
Science Operations Centre of the Solar Dynamic Observatory.

We use the 10.7 cm flux as a measure of solar activity (Tapping, 2013; Tapping
& Morton, 2013). These data are also public. 111Data are available from
http://www.spaceweather.gc.ca/solarflux/sx-en.php,
and https://omniweb.gsfc.nasa.gov/form/dx1.html.

For the purpose of looking at time variations and solar-activity related
changes, we use GONG data as one set and combine the MDI and HMI data to form
another set — if both sets show the same trends, we can be confident that we
are not seeing an instrumental effect. Since we are doing a differential
study, we choose one reference frequency set each for the ground-based and
space-based data. In each case, we choose sets from the minimum between Cycle
23 and Cycle 24. For the GONG sets, we use the averaged frequencies of GM 133,
136 and 139 as the reference set. These are non-overlapping sets and have an
average 10.7 cm flux of 68.1 SFU2221 solar flux unit or SFU is $10^{4}$ Jansky
$i.e.,\\\ 10^{-22}$ W m-2 Hz-1. For the space-based set, we use a specially
created 720-day data set spanning 2007.07.27 to 2009.07.16 as the reference
set. This set has an average 10.7 cm flux of 69.3 SFU

##  3 A clue to cycle-dependent structure changes

If the solar-cycle related changes of frequencies were a result only of
changes localized to the very shallow, near-surface layers of the Sun, the
changes would be a function of frequency alone, at least at low and
intermediate degrees; any degree differences can be accounted for by the
differences in the inertia of the modes. As in the case of all surface terms,
it is usual to fit an ad hoc, slowly varying function of frequency to model
the solar-cycle related change.

However, the form of the surface term need not be completely arbitrary.
Libbrecht & Woodard (1990) showed that the shifts that they found could be
modelled as the cube of the frequency divided by the mode inertia. This was
given some theoretical justification by Goldreich et al. (1991) who showed
that a photospheric perturbation could produce such a term. Gough (1990)
proposed an additional term, one that depends on the inverse of the frequency,
again scaled by the mode inertia. The cubic term could arise from the presence
of some factor that affects sound speed but not density, such as magnetic flux
tubes. The inverse term would reflect factors such as a change in the pressure
scale height near the surface. Ball & Gizon (2014) first used this to describe
the frequency dependence of the surface term in low-degree solar and
asteroseismic data.

The Ball & Gizon (2014) parameterization is given by:

$\delta\nu=\frac{1}{E_{nl}}\left[a_{-1}\left(\frac{\nu}{\nu_{\rm
ac}}\right)^{-1}+a_{3}\left(\frac{\nu}{\nu_{\rm ac}}\right)^{3}\right],$ (1)

where, $E_{nl}$ is the inertia of a model of degree $l$ and order $n$, and
$\nu_{\rm ac}$ is the acoustic cut-off frequency. Howe et al. (2017) showed
that this can be applied to both low-degree and intermediate degree solar
modes. What had not been examined directly in that work is how the fits behave
as a function of time, i.e., how the reduced $\chi^{2}$ behaves as a function
of time. The reduced $\chi^{2}$ for the frequency differences between the
individual GONG sets and the GONG reference set, as well as the frequency
differences of the MDI and HMI sets with respect to the MDI reference set, are
shown in Fig. 1. What is clear immediately is that there is a solar-cycle-
dependent change in the reduced $\chi^{2}$ and that the $\chi^{2}$ is largest
at the solar maxima.

Figure 1:  The variation of the reduced $\chi^{2}$ obtained by fitting Eq. 1
to the frequency differences are plotted as a function of time. The grey curve
in the background is the 10.7 cm flux whose values can be read from the axis
on the right.  Figure 2:  The time-variation of the average residuals for
modes with lower turning points between $0.7125\pm 0.0125$ R⊙ obtained from
the fits to Eq. 1. Only GONG data are shown for the sake of clarity. The grey
curve in the background is the 10.7 cm flux whose values can be read from the
axis on the right.

The variation of the reduced $\chi^{2}$ is an indication that there is more to
the solar cycle-dependence of the frequencies differences than just a
contribution from the near-surface layers. In Fig. 2 we show how the average
residuals for modes with lower-turning points between $0.7125\pm 0.0125$R⊙vary
as a function of time. Since modes that have their lower turning point within
any radius range can have different frequencies — recall that the classical
lower turning point of a mode is $r_{t}=c(r_{t})\sqrt{l(l+1)}/\omega$, where
$c(r_{t}$) is the sound speed at $r_{t}$ and $\omega=2\pi\nu$ — a frequency-
dependent solar-cycle related change alone cannot explain a solar-cycle
related change in the frequency of modes with a given lower turning point.
Fig. 2 points to changes in the structure around $0.7125\pm 0.0125$R⊙, which
means we cannot rule out changes at other radii either. We look for changes in
the internal sound speed by inverting the frequency difference between the
different data sets and the reference sets.

##  4 Inversions: Methods and systematic errors

Inversions for solar structure start with the linearized form of the
oscillation equations obtained using the variation principle:

$\frac{\delta\nu_{i}}{\nu_{i}}=\int K^{i}_{c^{2},\rho}\frac{\delta
c^{2}}{c^{2}}{\rm d}r+\int K^{i}_{\rho,c^{2}}\frac{\delta\rho}{\rho}{\rm
d}r+\frac{F(\nu_{i}}{E_{i}},$ (2)

where $\delta\nu_{i}/\nu_{i}$ is the relative frequency difference of the
$i$th mode, $\rho$ is the density, $K^{i}_{c^{2},\rho}$ and
$K^{i}_{\rho,c^{2}}$ are known functions of a reference solar model,
$F_{(}\nu_{i})$ is a slowly varying function function of frequency to account
for the surface term, and $E_{i}$ is the inertia of the $i$th mode. The aim of
inversions is to determine $\delta c^{2}/c^{2}$ or $\delta\rho/\rho$ given a
set of $\delta\nu/\nu$.

There are two primary ways of inverting Eq. 2 that have different aims, these
are The Regularised Least Squares (RLS) method, and the method of Optimally
Localised Averages (OLA). RLS aims to find the $\delta c^{2}/c^{2}$ and
$\delta\rho/\rho$ profiles that give the best fit to the data (i.e., give the
smallest residuals) while keeping uncertainties small; the aim of OLA is not
to fit the data at all, but to find linear combinations of the frequency
differences in such a way that the corresponding combination of kernels, the
averaging or resolution kernel, provides a localised average of the unknown
function, again while keeping uncertainties small. Details of how these
techniques are implemented can be found in Basu (2016). OLA and RLS inversions
are complementary in nature (Sekii 1997) and inversions can be trusted if both
inversion techniques return the same results.

There are two variants of OLA, the original version, which in helioseismic
literature is referred to as multiplicative OLA (MOLA) because a weighting
function used in the construction of the averaging kernels are multiplied with
the mode kernels, and subtractive OLA (SOLA; Pijpers & Thompson, 1992), where
the difference between the averaging kernels and a target kernel. SOLA has the
advantage that the target kernels can be selected for specific inversions,
such as the derivatives of the function rather than the function itself
(Pijpers & Thompson, 1994). Details of how SOLA and MOLA implementations
differ can be found in Basu (2016).

We use SOLA of the results obtained in this paper since it is faster and
allows us to directly invert for derivatives. We use mode kernels from Model S
(Christensen-Dalsgaard et al., 1996); our tests have shown that the kernels
from other standard solar models give identical results. In Fig. 3 we show the
inversion results for MDI set 3232 with respect to the 2-yr MDI set obtained
with all three methods, MOLA, SOLA and RLS, to determine the region where we
can trust the inversion results. As can be seen, all three methods agree for
radii below about $0.88$R⊙. We have not shown the results deeper than 0.6R⊙
because the RLS results there are linear, implying that they are dominated by
the regularization term rather than the data, and hence not to be trusted;
they also deviate from both SOLA and MOLA results.

Figure 3:  Inversion results for the difference in frequencies between MDI set
3232 and the 2year MDI set obtained using three different inversion methods.
Note that the results for the different methods agree at radii below about
0.88 R⊙. Although not shown in the figure, the match between the results
deteriorates below about 0.45R⊙.  Figure 4:  A few averaging kernels (i.e.,
resolution kernels) for the SOLA inversions shown in Fig. 3. Note that they
have very little structure away from the peak.

Since there is a trade-off between the resolution of the inversions and the
uncertainty in the results, we sacrificed resolution in order to increase the
precision of the results. The averaging kernels at a few different positions
obtained for the SOLA inversions are shown in Fig. 4. As can be seen, the
results at most radii will be correlated.

A possible source of systematic error is the structure of the averaging kernel
far from the target radius. This is caused mainly because of the lack of high
degree modes in the data set. In particular, the worry is that since the
averaging kernels are not strictly zero close to the surface, if there is a
large solar-cycle related change in the surface, this would leak into the
results at other radii giving a false impression about the deeper layers.
Recall that the SOLA inversion results are in effect (Basu, 2016)

$\left(\frac{\delta c^{2}}{c^{2}}\right)_{\rm
inv}=\sum_{i}c_{i}\frac{\delta\nu_{i}}{\nu_{i}}=\sum_{i}c_{i}\int
K^{i}_{c^{2},\rho}\frac{\delta c^{2}}{c^{2}}{\rm d}r,$ (3)

assuming that $\sum_{i}c_{i}\int K^{i}_{\rho,c^{2}}{\delta\rho/\rho}{\rm d}r$
and $\sum_{i}c_{i}F(\nu_{i}/E_{i}$, are small; $c_{i}$ are the inversion
coefficients. The term $\sum_{i}c_{i}K^{i}_{c^{2},\rho}$ is nothing but
$\mathcal{K}$, the averaging kernel, and $\sum_{i}c_{i}K^{i}_{\rho,c^{2}}$ is
usually called the ‘cross term’ kernel. Since the inverted value of the sound
speed at a radius $r_{0}$ is in effect $\int\mathcal{K}(r_{0};r)\delta
c^{2}/c^{2}{\rm d}r$, any structure in the averaging kernel away from $r_{0}$
will contribute to the solution. As can be seen from Fig. 4, the averaging
kernels are not zero close to the surface.

Figure 5:  Possible systematic error in solution caused by near-surface
structure in the averaging kernels. Top: the assumed difference in $\delta
c^{2}/c^{2}$, about 10% at the surface, with the inset showing the difference
near the helium ionization zone. Bottom: The systematic error at different
radii. The red squares assume that only $\delta c^{2}/c^{2}$ at the surface
matters; the blue triangles assume that there is also a relative difference in
density that is identical both in magnitude and sign to the assumed $\delta
c^{2}/c^{2}$, while the green circles assume that the density differences have
identical magnitude, but are opposite in sign to that of the assumed $\delta
c^{2}/c^{2}$.

To try and quantify the effect of the surface structure of the averaging and
cross-term kernels on the results, we assume a large change in structure (10%)
at the surface (see the top panel of Fig. 5, and we also assume that there are
changes as deep as the second helium ionization zone since there are indirect
detections of change in structure in those layers (Watson & Basu, 2020). We
then determine the effect of this change on the solutions at different radii.
We show the results in the lower panel of Fig. 5. We show the results assuming
that there is a difference only in sound speed, a change in both sound speed
and density with both changing the same way, and in one case where the changes
in sound speed and density have the same magnitude but opposite sign. As can
be seen in the figure, the systematic error because of unresolved surface
changes of 10% is less than a part in $10^{5}$ for radii below 0.8R⊙ and rises
thereafter. This will need to be kept in mind while interpreting the results
at different radii. Clearly, for us to claim that we see changes in $c^{2}$,
they need to be greater than 1 part in $10^{5}$.

##  5 Results

###  5.1 Changes in sound speed

In Figs. 6 and 7 we show the inversion results obtained with space based data
(MDI+HMI) and GONG respectively at three different radii. The results are
plotted both as a function of time and of the 10.7 cm flux. In the case of
GONG data, the results of all sets are plotted against time, but only non-
overlapping sets are plotted against the radio flux. It can be seen very
clearly that there is a time variation in the result. Note that the results
are with respect to the minimum between Cycle 23 and 24, and the results show
that the sound speed is lower at high activity than at low activity,
particularly in the convection zone. The changes are large enough that they
cannot be explained as systematic error due to changes in the near-surface
layers. Visually, the GONG results look more significant, but that is probably
because the sets overlap in time. The figures appear to show that there is a
variation in sound speed even below the base of the convection zone, but as
can be seen from Fig. 4, the averaging kernel for 0.65R⊙ samples the
convection zone too.

Figure 6:  Inversion results for MDI and HMI data sets at three radii. Left:
Results plotted as a function of time, with the curve in the background
showing the 10.7 cm flux plotted on a scale from 0 to 250 (not shown). The
horizontal line in each panel is the weighted average of the results. Note
that while the results at 0.65R⊙ and 0.85R⊙ are independent (see Fig. 4), the
result at 0.75R⊙ has an overlap with both 0.65 and 0.85R⊙. Right: The same
results plotted as a function of the 10.7 cm flux. The horizontal line is
again the weighted average. The red line shows a weighted linear least-
squares-fit to the points with the 95% confidence level shown as the dotted
lines.  Figure 7:  The same as Fig. 6, but for GONG data. Note that while we
plot all sets against time, we only plot non-overlapping sets against the
radio flux to avoid exaggerating the correlation.

As expected from the time variation, the changes are well correlated with the
10.7 cm flux at all radii. The correlation coefficients and their statistical
significance are listed in Table 1. Note that the lower the $p$ value, the
more the significance of the slope; a $p$-value of 0.05 corresponds to a 95%
confidence level. The slopes of the relation between the changes and the 10.7
cm flux obtained for the space-based and ground-based data do not completely
agree — GONG results are steeper. This is most likely due to the fact that the
modes sets used are different for GONG and MDI+HMI, and also that for the GONG
sets we have used a reference set of data that spans slightly less than a
year, while for MDI+HMI we subtract out a 2-year set. What is important
though, is that both sets show that there is an activity-related change in
structure. We also note that the change seems to be larger within the
convection zone than below it, but as mentioned earlier the results below the
convection zone have a substantial contribution from those in the convection
zone.

Table 1: Correlations and statistical significance Quantity | Radius | Pearson | Spearman |  $t$ value |  $p$ value
---|---|---|---|---|---
|  (R⊙) | coeff. | coeff. |  |
| MDI+HMI
$\delta c^{2}/c^{2}$ | 0.650 | $-0.54$ | $-0.46$ | $-5.83$ | $<0.00001$
$\delta c^{2}/c^{2}$ | 0.750 | $-0.71$ | $-0.69$ | $-10.63$ | $<0.00001$
$\delta c^{2}/c^{2}$ | 0.850 | $-0.80$ | $-0.74$ | $-12.14$ | $<0.00001$
$\Delta(\delta c^{2}/c^{2})$ | 0.713 | $0.37$ | $0.33$ | $3.91$ | $0.00015$
$\frac{\rm d}{{\rm d}r}(\delta c^{2}/c^{2})$ | 0.613 | $-0.02$ | $-0.04$ | $-0.48$ | $0.64$
$\frac{\rm d}{{\rm d}r}(\delta c^{2}/c^{2})$ | 0.713 | $-0.37$ | $-0.35$ | $-4.21$ | $0.00006$
$\frac{\rm d}{{\rm d}r}(\delta c^{2}/c^{2})$ | 0.813 | $-0.28$ | $-0.18$ | $-2.06$ | $0.04$
$(\delta c^{2}/c^{2})_{\rm alt}$ | 0.750 | $-0.61$ | $-0.60$ | $-8.17$ | $<0.00001$
| GONG
$\delta c^{2}/c^{2}$ | 0.650 | $-0.58$ | $-0.72$ | $-6.43$ | $<0.00001$
$\delta c^{2}/c^{2}$ | 0.750 | $-0.88$ | $-0.88$ | $-17.29$ | $<0.00001$
$\delta c^{2}/c^{2}$ | 0.850 | $-0.85$ | $-0.83$ | $-13.74$ | $<0.00001$
$\Delta(\delta c^{2}/c^{2})$ | 0.713 | $0.40$ | $0.37$ | $3.60$ | $0.00027$
$\frac{\rm d}{{\rm d}r}(\delta c^{2}/c^{2})$ | 0.613 | $-0.24$ | $-0.23$ | $-2.81$ | $0.33$
$\frac{\rm d}{{\rm d}r}(\delta c^{2}/c^{2})$ | 0.713 | $-0.45$ | $-0.50$ | $-4.71$ | $<0.00001$
$\frac{\rm d}{{\rm d}r}(\delta c^{2}/c^{2})$ | 0.813 | $-0.14$ | $0.05$ | $-0.45$ | $0.016$
$(\delta c^{2}/c^{2})_{\rm alt}$ | 0.750 | $-0.88$ | $-0.89$ | $-18.01$ | $<0.00001$
|  |  |  |  |

###  5.2 Differences below and above the convection zone

Figure 8:  The averaging kernels for inversions to determine the difference
above and below the convection-zone base. The vertical line marks the position
of the convection-zone base.  Figure 9:  Inversion results for the difference
in sound speed below and above the convection-zone base. As with previous
figures, although all GONG results are plotted as a function of time, only
results from non-overlapping data sets are plotted again the 10.7 cm flux.

The results in Figs. 6 and 7point to differences below and above the
convection zone, however, as mentioned earlier, the difference below the
convection zone are affected by differences above, and hence merely
subtracting the results will not give us the true difference. We make use of
the fact that for SOLA inversions we can define a suitable target averaging
kernel that will give us the quantity we seek, in this case, the difference in
the results at two locations. The averaging kernel obtained for directly
inverting for the difference in sound speed below and above the base of the
convection zone is shown in Fig. 8. Note that we get a clean averaging kernel
with no overlap below and above the base of the convection zone. The results
of the inversion are shown in Fig. 9. There is indeed a small, time-varying
difference between the layer below and above the convection zone, with the
sound speed in the convection zone being lower than that below it when the Sun
is active. The differences are small; however, the systematic error in the
results for an artificial profile like that in Fig. 5 is at most 3 parts in
$10^{6}$, i.e., smaller than the changes seen.

###  5.3 Is there a change in the position of the convection-zone base?

Figure 10:  (a) The relative squared sound-speed difference between models
with different positions of the base of the convection zone and a model with
the convection-zone base at 0.713R⊙. (b) The radial gradient of the $\delta
c^{2}/c^{2}$ shown in panel (a). (c) Results that would be obtained if we were
to invert the frequency differences between the models. These results are
obtained by multiplying the averaging kernel for each of the points with the
profiles in panel (b). Note that the value is highest at the base of the
convection zone and quickly tapers out on each side.  Figure 11:  The
averaging kernels for inversions to determine the derivative of $\delta
c^{2}/c^{2}$ at 0.613, 0.713 and 0.813R⊙.

The base of the convection zone is marked by a change in the temperature
gradient, consequently. the sound-speed difference for two models with
different $r_{b}$, the position of the convection-zone base shows a sharp rise
and a steep gradient near $r_{b}$, as is shown in Fig. 10.

We use the SOLA technique to determine the derivative of $\delta c^{2}/c^{2}$
directly; these are shown for three different radii in Fig. 11. Since the
averaging kernels are quite wide, we cannot expect to see the steep gradients
that are seen in Fig. 10(b), instead, the derived results will be smaller and
smoother, as shown in Fig. 10(c). However, Fig. 10(c) shows that even though
the inversions will smooth out the underlying gradient, the obtained
difference will be largest near the position of the convection-zone base.

Figure 12:  The radial derivative of $\delta c^{2}/c^{2}$ at 0.713R⊙ plotted
as a function of time and the 10.7 cm flux.

The inversion results at 0.713R⊙(the base of the solar convection zone as per
Christensen-Dalsgaard et al. 1991, Basu 1998) can be seen in Fig. 12. The
results appear to show a variation with time and activity, however, the
results are noisy enough that to actually determine the amount of change in
the position of the convection-zone base is extremely difficult. We do not
show the results at other radii since they are not very statistically
significant. Additionally, as seen from Fig. 10, if there is a change in the
position of the convection-zone base, it will be the largest in the region
around it. If we try to interpret Fig. 12(c) at face value with the aid of
Fig. 10, it would appear that the convection zone becomes deeper at activity
maximum.

###  5.4 Is the surface term fooling us?

Figure 13:  Inversion results at 0.750R⊙ obtained using an alternative surface
term. The slopes are almost identical to those obtained with the BG14 surface
term (Figs. 6 and 7).

Unlike earlier attempts in determining solar cycle-related structural changes,
we have used a very specific form of the surface term, that given in Eq. 1.
Thus one could wonder whether the results are simply an artefact of using that
form of the surface term, and perhaps that form does not really account for
all near-surface effects. To examine whether this is indeed the case, we
repeated the inversions with a conventional surface term. The results at
0.75R⊙ are shown in Fig. 13, and as can be seen, we see the same type of time
variation and the slope of the change against the 10.7 cm radio flux is almost
identical. Thus the form of the surface term is not causing the time-variation
that we find.

##  6 Discussion

We have used two solar cycle’s worth of helioseismic data to show that there
is indeed a solar-cycle related change in structure, particularly in the
convection zone. Watson & Basu (2020) had to use indirect means, we observe
the changes more directly by using conventional inversions, however, using
current mode sets, we can only probe changes between about 0.6 and 0.88 R⊙.
Although we started the study using the Ball & Gizon (2014) formulation of the
surface term, we find that the time dependence remains when we use the
convection forms of the surface term. The changes, although small, are
statistically significant.

Baldner & Basu (2008) had found a change $\delta c^{2}/c^{2}$ of $(7.23\pm
2.08)\times 10^{-5}$ near the base of the convection zone between the solar
minimum and maximum. The changes they find are localized at
$0.712^{+0.0097}_{-0.0029}$ R⊙. In contrast we get a change of $(2.56\pm
0.71)\times 10^{-5}$ using MDI data and $(3.99\pm 0.65)\times 10^{-5}$ using
GONG data. The discrepancy in the results is most likely due to differences in
the radial resolution of the two sets of inversions; in this work, we
sacrificed resolution to lower the uncertainties in the solution our
uncertainties are smaller nearly a factor of 3 compared with that of Baldner &
Basu (2008). Unlike the highly oscillatory nature of the Baldner & Basu (2008)
solution, our differences are smooth at all radii (see Fig. 3) which again
points to the effects of the radial resolution of the inversions.

The changes in solar dynamics between Cycles 23 and 24 were quite different
(Howe et al., 2017; Basu & Antia, 2019). In particular, changes in the
tachocline did not show the same activity-related change in the two cycles —
the change in the rotation rate across the tachocline was different for the
same value of the 10.7 cm flux. Given the uncertainty in our results, we are
not able to examine cycle-to-cycle changes between the response of internal
sound speed to solar activity. A principal component analysis of the two
cycles, done separately, might be able to detect cycle-to-cycle changes since
there is evidence that the frequency response to solar activity was different
in the two cycles (Broomhall & Nakariakov, 2015).

Our results show that as solar activity increases, the sound speed in the
region above the base of the convection zone, i.e., the tachocline region,
decreases compared with that below the convection-zone base. This result is
consistent with the assumption that the sound-speed changes are a result of
magnetic fields.

By directly inverting for the gradient of the sound-speed differences, we have
shown that there is a change in the gradient with the solar cycle. The change
is most significant at the base of the convection zone. Interpreted most
simply, this implies a change in the position of the convection-zone base.
However, with currently available frequencies, finding the change at any given
epoch is difficult. There had been earlier attempts to determine if the
position of the convection-zone base showed any change without much success
(Basu & Antia, 1999), however, those used a very limited data set (for
instance, Basu & Antia 1999 only had access to data from 1995 – 1998), and did
not use a differential method.

The results are intriguing enough to examine the results further. Given that
the current data sets are not ideal, the analyses need to be repeated with
frequencies obtained from time series that are longer; while longer time
series will smooth out some of the solar-cycle related changes, overlapping
sets will still reveal the changes. Two-year data sets at the solar maxima and
two years at the minima and a few two-year sets in between will help.

This work utilizes data from the National Solar Observatory Integrated
Synoptic Program, which is operated by the Association of Universities for
Research in Astronomy, under a cooperative agreement with the National Science
Foundation and with additional financial support from the National Oceanic and
Atmospheric Administration, the National Aeronautics and Space Administration,
and the United States Air Force. The GONG network of instruments is hosted by
the Big Bear Solar Observatory, High Altitude Observatory, Learmonth Solar
Observatory, Udaipur Solar Observatory, Instituto de Astrofísica de Canarias,
and Cerro Tololo Interamerican Observatory. This work also used data from the
Michelson Doppler Imager on board SOHO. SOHO is a project of international
cooperation between ESA and NASA. We also use data from the Helioseismic and
Magnetic Imager on board the Solar Dynamics Observatory and acknowledge the
HMI consortium for their support. Particular than are due to Timothy Larson
and Charles Baldner for analyzing the 2-year MDI data set. We would also like
to thank the referee for pointing out issues with the text. GONG, MDI, HMI,
Penticton Radio Observatory

## References

  * Antia & Basu (2013) Antia, H. M., & Basu, S. 2013, in Journal of Physics Conference Series, Vol. 440, Journal of Physics Conference Series, 012018
  * Baldner & Basu (2008) Baldner, C. S., & Basu, S. 2008, ApJ, 686, 1349, doi: 10.1086/591514
  * Ball & Gizon (2014) Ball, W. H., & Gizon, L. 2014, A&A, 568, A123, doi: 10.1051/0004-6361/201424325
  * Basu (1997) Basu, S. 1997, MNRAS, 288, 572
  * Basu (1998) —. 1998, MNRAS, 298, 719, doi: 10.1046/j.1365-8711.1998.01690.x
  * Basu (2002) Basu, S. 2002, in ESA Special Publication, Vol. 508, From Solar Min to Max: Half a Solar Cycle with SOHO, ed. A. Wilson, 7–14.  https://arxiv.org/abs/astro-ph/0205049
  * Basu (2016) —. 2016, Living Reviews in Solar Physics, 13, 2, doi: 10.1007/s41116-016-0003-4
  * Basu & Antia (1999) Basu, S., & Antia, H. M. 1999, in SOHO-9 Workshop on Helioseismic Diagnostics of Solar Convection and Activity, Vol. 9, 38
  * Basu & Antia (2000) Basu, S., & Antia, H. M. 2000, Solar Phys., 192, 449
  * Basu & Antia (2019) Basu, S., & Antia, H. M. 2019, ApJ, 883, 93, doi: 10.3847/1538-4357/ab3b57
  * Basu et al. (2007) Basu, S., Antia, H. M., & Bogart, R. S. 2007, ApJ, 654, 1146, doi: 10.1086/509251
  * Basu & Mandel (2004) Basu, S., & Mandel, A. 2004, ApJ, 617, L155, doi: 10.1086/427435
  * Broomhall & Nakariakov (2015) Broomhall, A. M., & Nakariakov, V. M. 2015, Sol. Phys., 290, 3095, doi: 10.1007/s11207-015-0728-6
  * Charbonneau et al. (1999) Charbonneau, P., Christensen-Dalsgaard, J., Henning, R., et al. 1999, ApJ, 527, 445
  * Chou & Serebryanskiy (2005) Chou, D.-Y., & Serebryanskiy, A. 2005, ApJ, 624, 420, doi: 10.1086/428925
  * Christensen-Dalsgaard et al. (1991) Christensen-Dalsgaard, J., Gough, D. O., & Thompson, M. J. 1991, ApJ, 378, 413, doi: 10.1086/170441
  * Christensen-Dalsgaard et al. (1996) Christensen-Dalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286, doi: 10.1126/science.272.5266.1286
  * Dziembowski et al. (2000) Dziembowski, W. A., Goode, P. R., Kosovichev, A. G., & Schou, J. 2000, ApJ, 537, 1026, doi: 10.1086/309056
  * Eff-Darwich et al. (2002) Eff-Darwich, A., Korzennik, S. G., Jiménez-Reyes, S. J., & Pérez Hernández, F. 2002, ApJ, 580, 574, doi: 10.1086/343063
  * Elsworth et al. (1990) Elsworth, Y., Howe, R., Isaak, G. R., McLeod, C. P., & New, R. 1990, Nature, 345, 322
  * Evans & Roberts (1992) Evans, D. J., & Roberts, B. 1992, Nature, 355, 230, doi: 10.1038/355230a0
  * Goldreich et al. (1991) Goldreich, P., Murray, N., Willette, G., & Kumar, P. 1991, ApJ, 370, 752, doi: 10.1086/169858
  * Gough (1990) Gough, D. O. 1990, Comments on Helioseismic Inference, ed. Y. Osaki & H. Shibahashi, Vol. 367, 283, doi: 10.1007/3-540-53091-6
  * Hill et al. (1996) Hill, F., Stark, P. B., Stebbins, R. T., et al. 1996, Science, 272, 1292
  * Howe (2009) Howe, R. 2009, Living Reviews in Solar Physics, 6, 1, doi: 10.12942/lrsp-2009-1
  * Howe et al. (2017) Howe, R., Basu, S., Davies, G. R., et al. 2017, MNRAS, 464, 4777, doi: 10.1093/mnras/stw2668
  * Howe et al. (2018a) Howe, R., Chaplin, W. J., Davies, G. R., et al. 2018a, MNRAS, 480, L79, doi: 10.1093/mnrasl/sly124
  * Howe et al. (2013) Howe, R., Christensen-Dalsgaard, J., Hill, F., et al. 2013, in Astronomical Society of the Pacific Conference Series, Vol. 478, Fifty Years of Seismology of the Sun and Stars, ed. K. Jain, S. C. Tripathy, F. Hill, J. W. Leibacher, & A. A. Pevtsov, 303
  * Howe et al. (2018b) Howe, R., Hill, F., Komm, R., et al. 2018b, ApJ, 862, L5.  https://arxiv.org/abs/1807.02398
  * Howe et al. (1999) Howe, R., Komm, R., & Hill, F. 1999, ApJ, 524, 1084
  * Howe et al. (2002) Howe, R., Komm, R. W., & Hill, F. 2002, ApJ, 580, 1172, doi: 10.1086/343892
  * Jain et al. (2009) Jain, K., Tripathy, S. C., & Hill, F. 2009, ApJ, 695, 1567, doi: 10.1088/0004-637X/695/2/1567
  * Kosovichev (1996) Kosovichev, A. G. 1996, ApJ, 469, L61
  * Kosovichev & Pipin (2019) Kosovichev, A. G., & Pipin, V. V. 2019, ApJ, 871, L20
  * Libbrecht & Woodard (1990) Libbrecht, K. G., & Woodard, M. F. 1990, Nature, 345, 779, doi: 10.1038/345779a0
  * Nishizawa & Shibahashi (1995) Nishizawa, Y., & Shibahashi, H. 1995, in Astronomical Society of the Pacific Conference Series, Vol. 76, GONG 1994. Helio- and Astro-Seismology from the Earth and Space, ed. R. K. Ulrich, J. Rhodes, E. J., & W. Dappen, 280
  * Palle et al. (1989) Palle, P. L., Regulo, C., & Roca Cortes, T. 1989, A&A, 224, 253
  * Pijpers & Thompson (1992) Pijpers, F. P., & Thompson, M. J. 1992, A&A, 262, L33
  * Pijpers & Thompson (1994) —. 1994, A&A, 281, 231
  * Scherrer et al. (1995) Scherrer, P. H., Bogart, R. S., Bush, R. I., et al. 1995, Sol. Phys., 162, 129
  * Scherrer et al. (2012) Scherrer, P. H., Schou, J., Bush, R. I., et al. 2012, Sol. Phys., 275, 207
  * Tapping (2013) Tapping, K. F. 2013, Space Weather, 11, 394
  * Tapping & Morton (2013) Tapping, K. F., & Morton, D. C. 2013, in Journal of Physics Conference Series, Vol. 440, Journal of Physics Conference Series, 012039
  * Watson & Basu (2020) Watson, C. B., & Basu, S. 2020, ApJ, 903, L29, doi: 10.3847/2041-8213/abc348
  * Woodard et al. (1991) Woodard, M. F., Kuhn, J. R., Murray, N., & Libbrecht, K. G. 1991, ApJ, 373, L81, doi: 10.1086/186056
  * Woodard & Noyes (1985) Woodard, M. F., & Noyes, R. W. 1985, Nature, 318, 449


A_{3}&A_{4}\end{matrix}\right)\left(\begin{matrix}\Lambda&0\\\
0&0\end{matrix}\right)V{\bf\omega}=\left(\begin{matrix}\Lambda&0\\\
0&0\end{matrix}\right)VB{\bf\omega}=\alpha\left(\begin{matrix}\Lambda&0\\\
0&0\end{matrix}\right)V{\bf\omega},$

since $B{\bf\omega}=\alpha{\bf\omega}$. Let ${\bf
v}\in\operatorname{M}_{r\times 1}(L)$ denote the first $r$ entries of
$V{\bf\omega}\in\operatorname{M}_{k\times 1}(L)$ and let ${\bf w}=\Lambda{\bf
v}$. Note that ${\bf v}$ is non-zero, since ${\bf\omega}$ is a basis and $V$
is non-singular. Also, ${\bf w}$ is non-zero, since
$\Lambda=\operatorname{diag}(\lambda_{1},\ldots,\lambda_{r})$ is non-singular.
Then (9.2) implies

$A_{1}{\bf w}=\alpha{\bf w}.$

Since ${\bf w}$ is non-zero, $\alpha$ is an eigenvalue of $A_{1}$. Hence,
$k=r$, $h_{1}$ is the characteristic polynomial of $A_{1}$, and $h_{2}$ is the
characteristic polynomial of $A_{4}$. Applying the induction process on $n$,
the statement of the theorem holds for $A_{4}\in\operatorname{M}_{n-k}(R)$ and
therefore, holds for $A$. ∎

## References

  * [BLP19] J. Bergfalk, M. Lupini, A. Panagiotopoulos, Definable $($co$)$homology, pro-torus rigidity, and $($co$)$homological classification, September 2019.
  * [CHL13] A. Clark, S. Hurder, and O. Lukina, Classifying matchbox manifolds, Geom. Topol. 23 (2019), no. 1, 1–27.
  * [D37] D. Derry, Uber Eine Klasse von Abelschen Gruppen, Proc. London Math. Soc. (2) 43 (1937), no. 7, 490–506.
  * [F73] L. Fuchs, Infinite abelian groups. Vol. II. Pure and Applied Mathematics. Vol. 36-II. Academic Press, New York-London, 1973.
  * [GPS19] T. Giordano, I. Putnam, C. Skau, ${\mathbb{Z}}^{d}$-odometers and cohomology, Groups Geom. Dyn. 13 (2019), no. 3, 909–938.
  * [GM81] R. Gütling, O. Mutzbauer, Classification of torsion free abelian minimax groups, Rend. Sem. Mat. Univ. Padova 64 (1981), 193–199.
  * [LMFDB] The LMFDB Collaboration, The L-functions and modular forms database, home page of the number field 4.0.144400.1, https://www.lmfdb.org/NumberField/4.0.144400.1, 2022, [Online; accessed 8 June 2022].
  * [M65] M. C. McCord, Inverse limit sequences with covering maps, Trans. Amer. Math. Soc. 114, 1965, 197–209.
  * [N72] M. Newman, Integral matrices, Pure and Applied Mathematics, Vol. 45. Academic Press, New York-London, 1972.
  * [SAGE] William A. Stein et al. Sage Mathematics Software (Version 9.6), The Sage Development Team, 2022, http://www.sagemath.org.
  * [S22] M. Sabitova, Generalized ideal classes in application to toroidal solenoids, Pacific Journal of Mathematics 318-1 (2022), 189–228.
  * [T49] O. Taussky, On a theorem of Latimer and MacDuffee, Canad. J. Math. 1 (1949), 300–302.
  * [V27] L. Vietoris, Über den höheren Zusammenhang kompakter Räume und eine Klasse von zusammenhangstreuen Abbildungen, Math. Ann. 97 (1927), no. 1, 454–472.