# On the inconsistency of $\ell_{1}$-penalised sparse precision matrix
estimation

Otte Heinävaara
Helsinki Institute for Information Technology HIIT
Department of Computer Science
University of Helsinki
&Janne Leppä-aho
Helsinki Institute for Information Technology HIIT
Department of Computer Science
University of Helsinki
Jukka Corander
Helsinki Institute for Information Technology HIIT
Department of Mathematics and Statistics
University of Helsinki
&Antti Honkela
Helsinki Institute for Information Technology HIIT
Department of Computer Science
University of Helsinki

###### Abstract

Various $\ell_{1}$-penalised estimation methods such as graphical lasso and
CLIME are widely used for sparse precision matrix estimation. Many of these
methods have been shown to be consistent under various quantitative
assumptions about the underlying true covariance matrix. Intuitively, these
conditions are related to situations where the penalty term will dominate the
optimisation. In this paper, we explore the consistency of $\ell_{1}$-based
methods for a class of sparse latent variable -like models, which are strongly
motivated by several types of applications. We show that all $\ell_{1}$-based
methods fail dramatically for models with nearly linear dependencies between
the variables. We also study the consistency on models derived from real gene
expression data and note that the assumptions needed for consistency never
hold even for modest sized gene networks and $\ell_{1}$-based methods also
become unreliable in practice for larger networks.

##  1 INTRODUCTION

Estimating the sparse precision matrix, i.e. the inverse covariance matrix,
from data is a very widely used method for exploring the dependence structure
of continuous variables. The motivation for the approach stems from the fact
that for a Gaussian Markov random field model, zeros in the precision matrix
translate exactly to absent edges in the corresponding undirected Gaussian
graphical model, thus being informative about the marginal and conditional
independence relationships among the variables.

The full $p$-dimensional covariance matrix contains $p(p+1)/2$ parameters,
making its accurate estimation from limited data difficult. Additionally, the
structure learning requires the inverse of the covariance, and matrix
inversion is in general a very fragile operation. To make the problem
tractable, some form of regularisation is typically needed. Direct
optimisation of the sparse structure would easily lead to very difficult
combinatorial optimisation problems. To avoid these computational
difficulties, several convex $\ell_{1}$-penalty-based approaches have been
proposed. Popular examples include $\ell_{1}$-penalised maximum likelihood
estimation (Meinshausen and Bühlmann, 2006), which also forms the basis for
the highly popular graphical lasso (glasso) algorithm (Friedman et al., 2008).
$\ell_{1}$ regularisation has also been used for example in a non-
probabilistic alternative with linear-programming-based constrained $\ell_{1}$
minimisation (CLIME) algorithm of Cai et al. (2011).

At the heart of the optimisation problems considered by all these methods is a
term depending on the $\ell_{1}$ norm of the estimated precision matrix.
$\ell_{1}$-penalisation-based approaches such as lasso are popular for sparse
regression, but they have a known weakness: in addition to promoting sparsity
they also push true non-zero elements toward zero (Zhao and Yu, 2006). In the
context of precision matrix estimation this effect would be expected to be
especially strong when some elements of the precision matrix are large, which
happens for scaled covariance matrices when the covariance matrix becomes ill-
conditioned. This phenomenon occurs frequently under the circumstances where
some of the variables are nearly linearly dependent.

In this paper we demonstrate a drastic failure of the $\ell_{1}$ penalised
sparse covariance estimation methods for a class of models that have a linear
latent variable structure where some variables depend linearly on others. For
such models even in the limit of infinite data, popular $\ell_{1}$ penalised
methods cannot yield results that are significantly better than based on
random guessing on any setting of the regularisation parameter. Yet these
models have a very clear sparse structure that becomes obvious from the
empirical precision matrix with an increasing $n$.

Given the huge popularity and success of linear models in modelling data,
structures like the one considered in our work are natural for various real
world data sets. Motivated by our discovery, we also explore the inconsistency
of $\ell_{1}$ penalised methods on models derived from real gene expression
data and find them poorly suited for such applications.

##  2 STRUCTURE LEARNING OF GAUSSIAN GRAPHICAL MODELS

###  2.1 BACKGROUND

We start with a quick recap on the basics of Gaussian graphical models in
order to formulate the problem of structure learning. For a more comprehensive
treatment of the subject, we refer to (Whittaker 1990; Lauritzen 1996). Let
$\mathbf{X}=(X_{1},\ldots,X_{p})^{\prime}$ denote a random vector following a
multivariate normal distribution with zero mean and a covariance matrix
$\mathbf{\Sigma}$, $\mathbf{X}\sim N_{p}(\mathbf{0},\mathbf{\Sigma})$. Let
$G=(V,E)$ be an undirected graph, where the $V=\\{1,\ldots,p\\}$ is the set of
nodes and $E\subset V\times V$ stands for the set of edges. The nodes in the
graph represent the random variables in the vector $\mathbf{X}$ and absences
of the edges in the graph correspond conditional independence assertions
between these variables. More in detail, we have that $(i,j)\not\in E$ and
$(j,i)\not\in E$ if and only if $X_{i}$ is conditionally independent of
$X_{j}$ given the remaining variables in $\mathbf{X}$.

In the multivariate normal setting, there is a one-to-one correspondence
between the missing edges in the graph and the off-diagonal zeros of the
precision matrix $\mathbf{\Omega}=\mathbf{\Sigma}^{-1}$, that is,
$\omega_{ij}=0\Leftrightarrow X_{i}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}X_{j}\ |\mathbf{X}\setminus\\{X_{i},X_{j}\\}$
(see, for instance, Lauritzen 1996, p. 129). Given an undirected graph $G$, a
Gaussian graphical model is defined as the collection of multivariate normal
distributions for $\mathbf{X}$ satisfying the conditional independence
assertions implied by the graph $G$.

Assume we have a complete (no missing observations) i.i.d. sample
$\mathbf{x}=(\mathbf{x}_{1},\ldots,\mathbf{x_{n}})$ from the distribution
$N_{p}(\mathbf{0},\mathbf{\Sigma})$. Based on the sample $\mathbf{x}$, our
goal in structure learning is to find the graph $G$, or equivalently, learn
the zero-pattern of $\mathbf{\Omega}$. The usual assumption is that the
underlying graph is sparse. A naive estimate for $\mathbf{\Omega}$ by
inverting the sample covariance matrix is practically never truly sparse for
any real data. Furthermore, if $n<p$ the sample covariance matrix is rank-
deficient and thus not even invertible.

One common approach to overcome these problems is to impose an additional
$\ell_{1}$-penalty on the elements of $\mathbf{\Omega}$ when estimating it.
This kind of regularisation effectively forces some of the elements of
$\mathbf{\Omega}$ to zero, thus resulting in sparse solutions. In the context
of regression models, this method applied on the regression coefficients goes
by the name of lasso (Tibshirani, 1996). There exists a wide variety of
methods making use of $\ell_{1}$-regularisation in the setting of Gaussian
graphical model structure learning (Yuan and Lin 2007; Meinshausen and
Bühlmann 2006; Banerjee et al. 2008; Friedman et al. 2008; Peng et al. 2009;
Cai et al. 2011; Hsieh et al. 2014).

###  2.2 $\ell_{1}$-REGULARISED METHODS

In this section we provide a brief review of selected examples of different
types of $\ell_{1}$-penalised methods.

####  2.2.1 Glasso

We begin with the widely used graphical lasso-algorithm (glasso) by Friedman
et al. (2008). Glasso-method maximises an objective function consisting of the
Gaussian log-likelihood and an $\ell_{1}$-penalty:

$\log\det(\mathbf{\Omega})-\textnormal{trace}(\mathbf{\Omega
S})-\lambda||\mathbf{\Omega}||_{1},$ (1)

where $\mathbf{S}$ denotes the sample covariance matrix and $\lambda>0$ is the
regularisation parameter controlling the sparsity of the solution. The
$\ell_{1}$ penalty, $||\boldsymbol{\Omega}||_{1}=\sum_{i,j}|\omega_{ij}|$, is
applied on all the elements of $\boldsymbol{\Omega}$, but the variant where
the diagonal elements are omitted is also common. The objective function (1)
is maximised over all positive definite matrices $\mathbf{\Omega}$ and the
optimisation is carried out in practice using a block-wise coordinate descent.

####  2.2.2 CLIME

Cai et al. (2011) approach the problem of sparse precision matrix estimation
from a slightly different perspective. Their CLIME-method (Constrained
$\ell_{1}$-minimisation for Inverse Matrix Estimation) seeks matrices
$\boldsymbol{\Omega}$ with a minimal $\ell_{1}$-norm under the following
constraint

$|\mathbf{S}\boldsymbol{\Omega}-\mathbf{I}|_{\infty}\leq\lambda,$ (2)

where $\lambda$ is the tuning parameter and
$|\mathbf{A}|_{\infty}=\max_{i,j}|a_{ij}|$ is the element-wise maximum. The
optimisation problem $\min_{\boldsymbol{\Omega}}||\boldsymbol{\Omega}||_{1}$
subject to the constraint (2) does not explicitly force the solution to be
symmetric, which is resolved by picking from estimated values $\omega_{ij}$
and $\omega_{ji}$ the one with a smaller magnitude into the final solution. In
practice, the optimisation problem is decomposed over variables into $p$ sub-
problems which are then efficiently solved using linear programming.

####  2.2.3 SCIO

Liu and Luo (2015) introduced recently a method called Sparse Column-wise
Inverse Operator (SCIO). The SCIO-method decomposes the estimation of
$\boldsymbol{\Omega}$ into a following smaller problems

$\min_{\boldsymbol{\beta}_{i}\in\mathbb{R}^{p}}\left\\{\frac{1}{2}\boldsymbol{\beta}_{i}^{T}\mathbf{S}\boldsymbol{\beta}_{i}-\textbf{e}_{i}^{T}\boldsymbol{\beta}_{i}+\lambda||\boldsymbol{\beta}_{i}||_{1}\right\\},$

where $\mathbf{S}$ and $\lambda$ are defined as before and $\textbf{e}_{i}$ is
an $i$:th standard unit vector. The regularisation parameter $\lambda$ can in
general vary with $i$ but this is omitted in our notation. The solutions
$\hat{\boldsymbol{\beta}}_{i}$ form the columns for the estimate of
$\boldsymbol{\Omega}$. Also for SCIO, the symmetry of the resulting precision
matrix must be forced, and this is done as described in the case of CLIME.

###  2.3 ALTERNATIVE METHODS

####  2.3.1 The naive approach

In addition to the above-mentioned $\ell_{1}$-penalised methods, we consider
two alternative approaches. In a ”naive” approach, we simply take the sample
covariance matrix, invert it, and then threshold the resulting matrix to
obtain a sparse estimate for the precision matrix. The threshold value is
chosen using the ground truth graph so that the naive estimator will have as
many non-zero entries as there are edges in the true graph. Setting the
threshold value according to the ground truth is of course unrealistic,
however, it is nevertheless interesting to compare the accuracy of this simple
procedure to the performance of the more refined $\ell_{1}$-methods, when also
their tuning parameters are chosen in a similar fashion.

####  2.3.2 FMPL

Lastly, we consider a Bayesian approach which is based on finding a graph with
a highest fractional marginal pseudo-likelihood (FMPL) by Leppä-aho et al.
(2016). The fractional marginal pseudo-likelihood is an approximation of the
marginal likelihood and it has been shown to be a consistent scoring function
in the sense that the true graph maximises it as the sample size tends to
infinity, under the assumption that data are generated from a multivariate
normal distribution. The FMPL-score decomposes over variables and in practice,
the method identifies optimal Markov blankets for each of the variables, which
are then combined into a proper undirected graph using any of the three
different schemes commonly employed in graphical model learning: OR, AND and
HC.

###  2.4 MODEL SELECTION CONSISTENCY

The assumptions required for a consistent model selection with an
$\ell_{1}$-penalised Gaussian log-likelihood have been studied, for instance,
in Ravikumar et al. (2011). The authors provide a number of conditions in the
multivariate normal model that are sufficient for the recovery of the zero
pattern of the true precision matrix $\boldsymbol{\Omega}^{*}$ with a high
probability when the sample size is large. For our purposes, the most relevant
condition is the following:

######  Assumption 1.

There exists $\alpha\in(0,1],$ such that

$\gamma:=||\Gamma_{S^{C}S}{(\Gamma_{SS})}^{-1}||_{\infty}\leq 1-\alpha.$ (3)

Here $S\subset V\times V$ is a set defining the support of
$\boldsymbol{\Omega}^{*}$, that is, the non-zero elements of
$\boldsymbol{\Omega}^{*}$ (diagonal and the elements corresponding to the
edges in the graphical model) and $S^{C}$ refers to the complement of $S$ in
$V\times V$. The $\Gamma$ term is defined via Kronecker product $\otimes$ as
$\Gamma=(\boldsymbol{\Omega}^{*})^{-1}\otimes(\boldsymbol{\Omega}^{*})^{-1}\in\mathbb{R}^{p^{2}\times
p^{2}}$ and $\Gamma_{AB}$ refers to the specific rows and columns of $\Gamma$
indexed by $A\subset V\times V$ and $B\subset V\times V$, respectively. The
norm in the equation is defined as $||A||_{\infty}=\max_{j}\sum_{i}|a_{ij}|$.

The above result applies to glasso. However, a quite similar result was
presented for SCIO in Liu and Luo (2015):

######  Assumption 2.

There exists $\alpha\in(0,1),$ such that

$\max_{1\leq i\leq
p}||\boldsymbol{\Sigma}^{*}_{\textbf{s}_{i}^{C}\textbf{s}_{i}}{(\boldsymbol{\Sigma}^{*}_{\textbf{s}_{i}\textbf{s}_{i}})}^{-1}||_{\infty}\leq
1-\alpha.$

Here $\boldsymbol{\Sigma}^{*}=(\boldsymbol{\Omega}^{*})^{-1}$ and
$\textbf{s}_{i}=\\{j\in\\{1,\ldots,p\\}\ |\ (\boldsymbol{\Omega}^{*})_{ij}\neq
0\\}$. Assumption 2 under the multivariate normality guarantees that the
support of $\boldsymbol{\Omega}^{*}$ is recovered by SCIO with a high
probability as the sample size gets large.

##  3 LATENT VARIABLE LIKE MODELS INDUCE INCONSISTENCY WITH
$\ell_{1}$-PENALISATION

Methods for sparse precision matrix estimation generally depend on an
objective function (such as log-likelihood) and a penalty function or
regulariser, which in a Bayesian setting is usually represented by the prior.
The ideal penalty function for many problems would be the $\ell_{0}$ “norm”
counting the number of non-zero elements: $||x||_{0}=\\#\\{i|x_{i}\neq 0\\}$.
This $\ell_{0}$ function is not a proper norm, but it provides a very
intuitive notion of sparsity. The main problem with its use is computational:
using $\ell_{0}$-penalisation leads to very difficult non-convex combinatorial
optimisation problems. The most common approach to avoid the computational
challenges is to use $\ell_{1}$-penalisation as a convex relaxation of
$\ell_{0}$. As mentioned above this works well in many cases but it comes with
a price, since in addition to providing the sparsity, $\ell_{1}$ also
regularises large non-zero values. Depending on the problem, as we demonstrate
here, this effect can be substantial and may cause $\ell_{1}$-regularised
methods to return totally meaningless results.

Intuitively, $\ell_{1}$-regularised methods are expected to fail when some
elements of the true precision matrix become so large that their contribution
to the penalty completely overwhelms the other parts of the objective and the
penalty. One example where this happens is when some set of variables depends
linearly on another set of variables. In such situation the covariance matrix
can become ill-conditioned and the elements of its inverse, the precision
matrix, grow. One example of when this happens is models with a linear latent
variable structure.

Figure 1: Left: Graphical representation of a latent variable model as an
undirected graphical model for a case with somewhat sparse $\mathbf{A}$.
Right: The adjacency matrix of the graph showing the sparse pattern of non-
zero elements in the corresponding precision matrix.

Let us consider a model for
$\mathbf{x}\in\mathbb{R}^{d_{1}},\mathbf{y}\in\mathbb{R}^{d_{2}}$, where
$\mathbf{y}=\mathbf{A}\mathbf{x}+\epsilon$. The graphical structure of the
model and the corresponding precision matrix structure are illustrated in Fig.
1. Assuming
$\mathbf{x}\sim\mathcal{N}(0,\sigma_{x}^{2}I),\epsilon\sim\mathcal{N}(0,\sigma_{\epsilon}^{2}I)$,
the covariance of the concatenated vectors
$(\mathbf{x}^{T},\mathbf{y}^{T})^{T}$ is given by the block matrix

$\mathrm{Cov}((\mathbf{x}^{T},\mathbf{y}^{T})^{T})=\mathbf{C}=\sigma_{x}^{2}\begin{pmatrix}I&\mathbf{A}^{T}\\\
\mathbf{A}&\mathbf{A}\mathbf{A}^{T}+\sigma_{\epsilon}^{2}I\end{pmatrix}.$ (4)

The covariance matrix has an analytical block matrix inverse (Lu and Shiou,
2002)

$\mathbf{C}^{-1}=\sigma_{x}^{-2}\begin{pmatrix}I+\sigma_{\epsilon}^{-2}\mathbf{A}^{T}\mathbf{A}&-\sigma_{\epsilon}^{-2}\mathbf{A}^{T}\\\
-\sigma_{\epsilon}^{-2}\mathbf{A}&\sigma_{\epsilon}^{-2}I\end{pmatrix}.$ (5)

This precision matrix recapitulates the conditional independence result for
Gaussian Markov random fields: the lower right block is diagonal because the
variables in $\mathbf{y}$ are conditionally independent of each other given
$\mathbf{x}$. The matrix is clearly sparse, so we would intuitively assume
sparse precision matrix estimation methods should be able to recover it. The
non-zero elements do, however, depend on $\sigma_{\epsilon}^{-2}$ which can
make them very large if the noise $\sigma_{\epsilon}^{2}$ is small.

It is possible to evaluate and bound the different terms of Eq. (1) evaluated
at the ground truth for these models:

$\displaystyle\log\det(\mathbf{C}^{-1})$
$\displaystyle=-d_{2}\log\sigma_{\epsilon}^{2}$
$\displaystyle-\textnormal{trace}(\mathbf{C}\mathbf{C}^{-1})$
$\displaystyle=-(d_{1}+d_{2})$ $\displaystyle-\lambda||\mathbf{C}^{-1}||_{1}$
$\displaystyle<-\lambda\sigma_{\epsilon}^{-2}(d_{2}+2||\mathbf{A}||_{1}).$

The magnitude of the last penalty term clearly grows very quickly as
$\sigma_{\epsilon}^{2}$ decreases. Clearly the magnitude of the two first log-
likelihood terms grows much more slowly as they only depend on
$\log\sigma_{\epsilon}^{2}$. Thus the total value of Eq. (1) decreases without
bound as $\sigma_{\epsilon}^{2}$ decreases.

Forgetting the ground truth, it is easy to see that one can construct an
estimate $\mathbf{\Omega}$ for which the objective remains bounded. If we
assume all values of $\mathbf{C}$ to be $\leq 1$ (after normalisation),

$\textnormal{trace}(\mathbf{C}\mathbf{\Omega})\leq||\mathbf{\Omega}||_{1}.$

As the other terms only depend on $\mathbf{\Omega}$ it is easy to choose
$\mathbf{\Omega}$ so that they remain bounded. The estimate $\mathbf{\Omega}$
that yields these values will in many cases not have anything to do with
$\mathbf{C}^{-1}$, as seen in the experiments below.

##  4 EXPERIMENTS

We tested the performance of glasso, SCIO and CLIME as well as FMPL using the
model structure introduced in Sec. 3. The performance of the methods was
investigated by varying the noise variance $\sigma_{\epsilon}^{2}$, and the
sample size $n$. The model matrix $\mathbf{A}$ was created as a
$(d_{2},d_{1})$-array of independent normal random variables with mean $0$ and
variance $1$. The majority of the tests were run using input dimensionality
$d_{1}=2$, output dimensionality $d_{2}=10$ and noise variance
$\sigma_{\epsilon}^{2}=0.1^{2}$ but we also tested varying these settings. For
each individual choice of noise and sample size, $k=50$ different matrices
$\mathbf{A}$ were generated and the results were averaged.

Generating $n$ samples using model described, data were normalised and
analysed using the five different methods. We calibrated the methods in a way
that number of edges in the resulting graph would match the true number.
Similarly, we thresholded the naive method by taking inverse matrix directly
to output the correct number of edges. The FMPL method has no direct tuning
parameters so we used its OR mode results as such. Similar tuning is not
possible in a real problem where the true number of edges is now known. The
tuning represents the best possible results the methods could obtain with an
oracle that provides an optimal regularisation parameter.

We evaluated the results using the Hamming distance between the ground truth
and the inferred sparsity pattern, i.e. the number of incorrect edges and non-
edges which were treated symmetrically. For methods returning the correct
number of edges, this value is directly related to the precision $pr$ through

$d_{\text{Hamming}}=2(1-pr)N_{\text{true positives}}$

or conversely

$pr=1-\frac{d_{\text{Hamming}}}{2N_{\text{true positives}}}.$

We will nevertheless use the Hamming distance as it enables fair comparison
with FMPL that sometimes returns a different number of edges.

Figure 2: Performances of different methods on latent variable like model with
100 samples. (Lower values are better.) Figure 3: Performances of different
methods on latent variable like model with 1000 samples. (Lower values are
better.)

Figs. 2 and 3 show the Hamming distance obtained by the different methods as a
function of the noise level when using 100 and 1000 samples, respectively. The
results show that especially for low but also for high noise levels, the
$\ell_{1}$-based methods all perform very poorly with especially glasso and
CLIME performing very close to random guessing level for low noise levels
$\sigma_{\epsilon}\leq 0.1$. The naive inverse and FMPL work much better up to
moderate noise levels of $\sigma_{\epsilon}\approx 2$ after which the noise
starts to dominate the signal and the performance of all methods starts to
drop. SCIO is a little better than the other $\ell_{1}$-based methods but
clearly worse than FMPL and naive in the low noise regime.

Figure 4: Performances of different methods on latent variable like model with
varying output dimensionality. (Lower values are better.) Figure 5:
Performances of different methods on latent variable like model with varying
input dimensionality. (Lower values are better.)

Fig. 4 shows the results when changing the output dimensionality $d_{2}$ from
10. The results show that the performance of all $\ell_{1}$-based methods is
very poor across all $d_{2}$. Glasso performance is close to random guessing
level across the entire range considered, while CLIME is slightly better for
$d_{2}\geq 18$ and SCIO slightly better across the entire range. Both FMPL and
naive are significantly better than any of the $\ell_{1}$-based methods.

Fig. 5 shows the corresponding result when changing the input dimensionality
$d_{1}$. The results are now quite different as all methods are better than
random especially for larger values. SCIO still outperforms CLIME which
outperforms glasso. FMPL is really accurate for small $d_{1}$ but degrades for
larger $d_{1}$ while the naive method is the most accurate in almost all
cases.

Figure 6: Contributions of the different terms of the glasso objective (1) for
latent variable like model with 1000 samples. The green curves show the
contributions of the first two terms of Eq. (1) and the blue curves show the
contributions of the last penalty term. Solid lines show the result of the
glasso optimal solution while dashed lines show the result for the true
solution.

To further illustrate the behaviour of glasso on these examples, Fig. 6 shows
the contributions of the different parts of the glasso objective function (1)
as a function of the noise level both for the true solution (“truth”) as well
as the glasso solution. The results show that for low noise levels the penalty
incurred by the true solution becomes massive. The glasso solution has a much
lower log-likelihood (“logl”) than ground truth but this is amply compensated
by the significantly smaller penalty. As the noise increases, the penalty of
the true solution decreases and the glasso solution converges to similar
values.

###  4.1 NECESSITY OF ASSUMPTION 1

It can be checked that the norm $\gamma$ in Assumption 1 and Eq. (3) for
latent-variable-like models depends on the scale of $\mathbf{A}$. We took
advantage of this by creating examples with different values of $\gamma$ and
testing the precision of glasso using the true covariance which corresponds to
infinite data limit. The results of this experiment are shown in Fig. 7. The
results verify that glasso consistently yields perfect results when $\gamma<1$
which is a part of the sufficient conditions for consistency of glasso. As
$\gamma$ grows and the sufficient conditions are no longer satisfied, it is
clearly seen that the accuracy of glasso starts to deteriorate rapidly. This
suggests that the sufficient condition of Assumption 1 is in practice also
necessary to ensure consistence.

Figure 7: Precision of glasso on infinite data as a function of the norm
$\gamma$ of Assumption 1 and Eq. (3). Values to the left of the green vertical
line satisfy this condition while values to the right violate it. (Higher
values are better.)

##  5 INCONSISTENCY FOR MODELS OF REAL GENE EXPRESSION DATA

We tested how often the problems presented above appear in real data using the
“TCGA breast invasive carcinoma (BRCA) gene expression by RNAseq
(IlluminaHiSeq)” data set (Cancer Genome Atlas Network, 2012) downloaded from
https://genome-cancer.ucsc.edu/proj/site/hgHeatmap/. The data set contains
gene expression measurements for 20530 genes for $n=1215$ samples. After
removing genes with a constant expression across all samples there are
$p=20252$ genes remaining.

In order to test the methods we randomly sampled subsets of $d$ genes and
considered the correlation matrix $\mathbf{C}_{0}$ over that subset. We
generated sparse models with known ground truth by computing the corresponding
precision matrix $\boldsymbol{\Lambda}_{0}$ from the empirical correlation
matrix, setting elements with absolute values below chosen cutoff $\delta=0.1$
to 0 to obtain

$\boldsymbol{\Lambda}_{ij}=\begin{cases}(\boldsymbol{\Lambda}_{0})_{ij}\quad&\text{if
}|(\boldsymbol{\Lambda}_{0})_{ij}|>\delta\\\
0\quad&\text{otherwise}\end{cases}$ (6)

and the testing covariance matrix $\mathbf{C}=\boldsymbol{\Lambda}^{-1}$. The
cutoff lead to networks that were sparse with on average 60% zeros in the
precision matrix.

Figure 8: Testing the condition of Assumption 1 of Ravikumar et al. (2011) in
Eq. (3) on real gene expression data showing the fraction of random subsets of
$d$ genes that fulfil the requirement and various relaxations. The condition
(3) requires $\gamma<1$, but the figure shows results also for larger $\gamma$
cutoffs.

Fig. 8 shows the fraction of covariances derived from random subsets of $d$
genes that satisfy the Assumption 1 of Ravikumar et al. (2011) ($c=1$) as well
as the fraction of values below more relaxed bounds. The figure shows that the
assumption is reliably satisfied only for very small $d$ while for $d\geq 20$,
the assumption is essentially never satisfied. Based on the results of Fig. 7
it is likely that glasso results will degrade significantly by for $\gamma>10$
and beyond which are very common for large networks.

We further studied how accurately glasso can recover the graphical structures
when the data were generated using the precision matrices described above. We
used a similar thresholding with a cut-off value of $0.1$ in order to first
form sparse precision matrices for a random subset of genes with given
dimension. These matrices were then inverted to obtain covariance matrices. We
checked that the resulting matrices were positive definite and then used them
to sample multivariate normal data with zero mean with different sample sizes.

The obtained data sets were centred and scaled before computing the sample
covariance which was used as input to the glasso algorithm. The regularisation
parameter was chosen with the aid of the ground truth graph, so that the the
graph identified by glasso would contain as many edges as there were in the
real graph. Results are shown in Figure 9. The results show that glasso
performance decreases as the network size increases and is approaching that of
random guessing for the largest networks considered here.

Figure 9: Average precisions for glasso with different dimensions and sample
sizes of the real gene expression data, higher values are better. The
precision obtained by random guessing is also illustrated. Figure 10: Average
contributions of the different terms of the glasso objective function (1) on
real gene expression data over random subsets of $d$ genes. The values are
shown for the $\ell_{1}$ penalty term as well as the unnormalised log-
likelihood, divided by $d$ to make them comparable. Solid lines show the
values for glasso result while dashed lines show the result for ground truth.

Fig. 10 shows the contributions of different parts of the glasso objective
function (1) as a function of the number of genes $d$. The regularisation
parameter $\lambda$ of glasso was tuned to return a solution with the same
number of edges as in the true solution. We used the glasso implementation of
scikit-learn (Pedregosa et al., 2011), which ignores the diagonal terms of
$\boldsymbol{\Omega}$ when computing the penalty. The figure shows clearly how
the penalty term for the true solution increases superlinearly as a function
of $d$. (A linear increase would correspond to a horizontal line.) The result
is even more striking given that the optimal $\lambda$ decreases slightly as
$d$ increases. The penalty contribution for glasso solution increases much
more slowly. The excess loss in log-likelihood from glasso solution increases
as $d$ increases, but this is compensated by a larger saving in the penalty.
Together these suggest that glasso solutions are likely to remain further away
from ground truth as $d$ increases.

##  6 DISCUSSION

The class of latent variable like models presented in Sec. 3 is an interesting
example of models that have a very clear sparse structure, which all
$\ell_{1}$-penalisation-based methods seem unable to recover even in the limit
of infinite data. This class complements the previously considered examples of
models where glasso is inconsistent including the “two neighbouring triangles”
model of Meinshausen (2008) and the star graph of Ravikumar et al. (2011), the
latter of which can be seen as a simple special case of our example.

An important question arising from our investigation is how significant the
discovered limitation to inferring sparse covariance matrices is in practice,
i.e. how common are the latent variable like structures in real data sets.
Given the popularity and success of linear models in diverse applications it
seems plausible such structures could often exist in real data sets, either as
an intrinsic property or as a result of some human intervention, e.g. through
inclusion of partly redundant variables.

The gene expression data set is a natural example of an application where
graphical model structure learning has been considered. The original glasso
paper (Friedman et al., 2008) contained an example on learning gene networks,
although from proteomics data. Other authors (e.g. Ma et al., 2007) have
applied Gaussian graphical models and even glasso (e.g. Menéndez et al., 2010)
to gene network inference from expression data. Our experiments on the TCGA
gene expression data suggest that in such applications it is advisable to
consider the conditions for the consistency of $\ell_{1}$ penalised methods
very carefully when planning to apply those.

Previous publications presenting new methods for sparse precision matrix have
typically tested the method on synthetic examples where the true precision
matrix is specified to contain mostly small values. Specifying the precision
matrix provides a convenient way to generate test cases as the sparsity
pattern can be defined very naturally through it. At the same time, this
excludes any models that have an ill-conditioned covariance. As shown by our
example, such ill-conditioned covariances arise very naturally from model
structures that are plausible from the application perspective.

Ultimately, our results suggest that users of the numerous $\ell_{1}$
penalised methods should be much more careful about checking whether the
conditions of consistency for precision matrix estimation are likely to be
fulfilled in the application area of interest.

#### Acknowledgements

This work was supported by the Academy of Finland [259440 to A.H., 251170 to
J.C.] and the European Research Council [239784 to J.C.].

## References

  * Banerjee et al. (2008) O. Banerjee, L. El Ghaoui, and A. d’Aspremont.  Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data.  _Journal of Machine Learning Research_ , 9:485–516, June 2008.
  * Cai et al. (2011) T. Cai, W. Liu, and X. Luo.  A constrained $\ell_{1}$ minimization approach to sparse precision matrix estimation.  _Journal of the American Statistical Association_ , 106(494):594–607, Jun 2011.
  * Cancer Genome Atlas Network (2012) Cancer Genome Atlas Network.  Comprehensive molecular portraits of human breast tumours.  _Nature_ , 490(7418):61–70, Oct 2012.
  * Friedman et al. (2008) J. Friedman, T. Hastie, and R. Tibshirani.  Sparse inverse covariance estimation with the graphical lasso.  _Biostatistics_ , 9(3):432–441, Jul 2008.
  * Hsieh et al. (2014) C. Hsieh, M. A. Sustik, I. S. Dhillon, and P. D. Ravikumar.  QUIC: quadratic approximation for sparse inverse covariance estimation.  _Journal of Machine Learning Research_ , 15(1):2911–2947, 2014.
  * Lauritzen (1996) S. Lauritzen.  _Graphical Models_.  Clarendon Press, 1996.  ISBN 9780191591228.
  * Leppä-aho et al. (2016) J. Leppä-aho, J. Pensar, T. Roos, and J. Corander.  Learning Gaussian graphical models with fractional marginal pseudo-likelihood.  _arXiv:1602.07863_ , 2016.
  * Liu and Luo (2015) W. Liu and X. Luo.  Fast and adaptive sparse precision matrix estimation in high dimensions.  _Journal of Multivariate Analysis_ , 135:153 – 162, 2015\.
  * Lu and Shiou (2002) T.-T. Lu and S.-H. Shiou.  Inverses of $2\times 2$ block matrices.  _Computers & Mathematics with Applications_, 43(1-2):119–129, Jan 2002.
  * Ma et al. (2007) S. Ma, Q. Gong, and H. J. Bohnert.  An Arabidopsis gene network based on the graphical Gaussian model.  _Genome Res_ , 17(11):1614–1625, Nov 2007.
  * Meinshausen (2008) N. Meinshausen.  A note on the Lasso for Gaussian graphical model selection.  _Statistics & Probability Letters_, 78(7):880–884, May 2008.
  * Meinshausen and Bühlmann (2006) N. Meinshausen and P. Bühlmann.  High-dimensional graphs and variable selection with the lasso.  _The Annals of Statistics_ , 34(3):1436–1462, Jun 2006.
  * Menéndez et al. (2010) P. Menéndez, Y. A. I. Kourmpetis, C. J. F. ter Braak, and F. A. van Eeuwijk.  Gene regulatory networks from multifactorial perturbations using Graphical Lasso: application to the DREAM4 challenge.  _PLoS One_ , 5(12):e14147, 2010.
  * 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\.
  * Peng et al. (2009) J. Peng, P. Wang, N. Zhou, and J. Zhu.  Partial correlation estimation by joint sparse regression models.  _Journal of the American Statistical Association_ , 104(486):735–746, 2009.
  * Ravikumar et al. (2011) P. Ravikumar, M. J. Wainwright, G. Raskutti, and B. Yu.  High-dimensional covariance estimation by minimizing $\ell_{1}$-penalized log-determinant divergence.  _Electronic Journal of Statistics_ , 5:935–980, 2011.
  * Tibshirani (1996) R. Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society, Series B_ , 58:267–288, 1996.
  * Whittaker (1990) J. Whittaker.  _Graphical Models in Applied Multivariate Statistics_.  John Wiley & Sons, 1990.
  * Yuan and Lin (2007) M. Yuan and Y. Lin.  Model selection and estimation in the Gaussian graphical model.  _Biometrika_ , 94(1):19–35, 2007.
  * Zhao and Yu (2006) P. Zhao and B. Yu.  On model selection consistency of lasso.  _Journal of Machine Learning Research_ , 7:2541–2563, 2006\.


# Study of Spherical Array Target for Long-Baseline Neutrino Experiment

Quynh M. Nguyen  University of Minnesota, Minneapolis, MN, USA, 55455  Fermi
National Accelerator Laboratory (FNAL), Batavia, IL, USA,60510  Alberto
Marchionni  Fermi National Accelerator Laboratory (FNAL), Batavia, IL,
USA,60510

(August 8, 2014)

###### Abstract

Spherical Array Target was studied by implementing the geometry in LBNE’s Beam
Simulation source code g4lbne version v3r2p4 and Monte Carlo. To compare with
Nominal LBNE target, unoscillated Far Detector neutrino flux was produced
using different parameters: sphere diameter 17mm and 13mm, different
longitudinal positions, two interaction length and less, beam size R/3 and
1.7mm, and beam offset from 50 $\mu$m to 1mm. The 1.86 interaction length
(901mm), 17mm diameter target, with beam size 1.7mm gives higher $\nu_{\mu}$
flux up to 10% from 0-3.5 GeV, and suppresses flux up to 70% at energy higher
than 3.5 GeV.

##  I Introduction to LBNE and the Target

The Long-Baseline Neutrino Experiment (LBNE) collaboration plans a
comprehensive experiment that will study CP violation, neutrino mass hierarchy
and other physics beyond the Standard Modellbne  using a high intensity 1300km
baseline accelerator neutrino beam and an advanced 34kton liquid argon TPC as
the far detector.

The LBNE Beamline at Fermilab will be designed to provide a neutrino beam of
sufficient intensity and proper energy range to meet the goals of the
experiment. The design is a conventional, horn-focused beamline volume2 .
Proton beam from Fermilab Main Injector will be extracted and transport to a
target area where collisions generate a beam of charged particles mainly pions
and kaons. The charge particles will be focused by magnetic horns, lead into
decay-pipe tunnel where they decay to generate neutrino beam, aiming toward
the Far Detector.

The current nominal target is a segmented graphite target, such as one used by
NuMi. The spherical array target studied and designed by Rutherford Appleton
Laboratory (RAL) promises interesting mechanical properties but there was no
study on its ability to produce neutrino flux at the Far Detector. This paper
presents a study of the spherical array target (SAT) for the LBNE’s beam-
simulation group. Geometry of the new target was developed in LBNE Geant4’s
derivative, g4lbne version v3r2p4.

##  II Design of SAT and Study Motivation

Figure 1: ral Concept design of an array of spheres with triple helix spacers
and Perspex outter casing

As designed by RAL, the target is an array of spheres, containing in a helical
spacer and a Titanium support cylindrical tube (Fig. 1). The spheres are made
of Beryllium. Properties of graphite and beryllium are provide in Table. 1.
Except a higher modulus of elasticity, Beryllium generally offers better
mechanical properties as a target material than graphite such as higher
specific heat, thermal conductivity and yield strength.

The spherical design is to avoid sharp edges and to offer most homogeneous
stress field. Furthermore, spheres also have no structural weak point ral  and
provide more surface for cooling. Stress simulations by RAL have shown that,
an array of spheres of diameter as small as 13mm can keep the mechanical
stress below the design point at beam power 2.3MW. (Fig. 3).

Figure 2: Dimension of sphere targets (13mm and 17 mm) compared to NuMi
target. Helium can size used in this study is 36mm.

Thus, spheres diameter of 13mm and 17mm were chosen. The diameter of a simple
Helium can was inferred from T2K’s target t2k  to be approximately 36mm. The
can is made of Titanium 0.4mm thick.

Figure 3: ral Peak stress for a 2 sigma off centre beam. Study done by RAL | Graphite (POCO ZXF-5Q) | Beryllium (S-65C)
---|---|---
Apparent density (g cm-3) | 1.81 | 1.848111Particle Data Group
Compressive Strength (MPa) | 195 | 260
Tensile Strength (MPa) | 90 | 370
Yield Strength (MPa) | - | 240 222GoodFellow.com
Modulus of Elasticity (GPa) | 12.5 | 310
Thermal Conductivity (W m-1 K-1) | 70 | 200
Coeff. of Thermal Expansion ($\mu$m m-1 K-1) | 8.1 | 10.7
Specific Heat (J kg-1 K-1) | 710 | 1770
Nuclear interaction length $\lambda_{I}$ 333Particle Data Group(g cm-2) | 85.8 | 77.8
Table 1: volume2  Properties of graphite and beryllium at 20∘C, from the
manufacturers, except as noted

##  III Implementation

###  III.1 Target length calculation

The length of SAT was calculated such that on average, the beam traversing SAT
longitudinally has a path length of two interaction lengths. The path length
of proton beams depends on beam’s offset from target’s center, and the beam is
assumed to have a bi-dimensional Gaussian distribution center at target
longitudinal axis with sigma equal 1/3 of the sphere radius, which means
99.89% of the beam is within the cross-section (However, smaller beam size was
also studied). The effective interaction length can be expressed as:

$N\times\int_{-\sqrt{R^{2}-x^{2}}}^{\sqrt{R^{2}-x^{2}}}\int_{-R}^{R}\\!\left(\frac{1}{{2\pi}\sigma_{x}\sigma_{y}}e^{-\frac{x^{2}}{2\sigma_{x}^{2}}-\frac{y^{2}}{2\sigma_{y}^{2}}}\right)\times
2\sqrt{R^{2}-x^{2}-y^{2}}\,\mathrm{d}x\mathrm{d}y=2\times\text{interaction
length}$ (1)

where $2\sqrt{R^{2}-x^{2}-y^{2}}$ is the parallel path length of a beam
traversing through a sphere (see Fig. 4), N is the number of spheres in the
target. Here, both $\sigma_{x}$ and $\sigma_{y}$ equal to R/3. Beryllium
interaction length is 421mm, calculated by the ratio Nuclear Interaction
length/Density from Table. 1. Solving for number of spheres give 57 spheres
(rounded) of diameter 17mm (total length 969mm) and 74 spheres of diameter
13mm (total length 975mm). In the source code, a tolerance of 5 microns is
also added as allowed gap between each spheres.

Figure 4: Path length of off-centre beam. Off-centre beam passes through less
material. Thus the target is longer than a simple cylinder to make two
effective interaction lengths.

###  III.2 Probability of interaction of off-axis proton

When the beam is offset, it traverses less target material, so does the
probability of interaction. For a Gaussian beam distribution, the probability
of interaction for off-centre beam is:

$P(\mu_{x},\mu_{y})=\int_{-\sqrt{R^{2}-x^{2}}}^{\sqrt{R^{2}-x^{2}}}\int_{-R}^{R}\\!\left(\frac{1}{{2\pi}\sigma_{x}\sigma_{y}}e^{-\frac{(x-\mu_{x})^{2}}{2\sigma_{x}^{2}}-\frac{(y-\mu_{y})^{2}}{2\sigma_{y}^{2}}}\right)\times\left(1-e^{-\frac{2N\sqrt{R^{2}-x^{2}-y^{2}}}{\text{interaction
length}}}\right)\,\mathrm{dx}\mathrm{dy}$ (2)

In study of beam position using Monte Carlo simulations, beam offset was set
in x direction (y direction would give the same result regarding the
simulation). To check on the Monte Carlo, the integral in Eq. 2 with $\mu_{x}$
as beam offset, $\mu_{y}=0$, is used to approximate $P(\mu_{x},\mu_{y})$, in
comparison to on-centre beam. Relative magnitude of P $(\mu_{x},\mu_{y})$
should agree with neutrino flux.

Beam offset (mm) | 0 | 0.05 | 0.1 | 0.2 | 0.4 | 1 | 2 | 3
---|---|---|---|---|---|---|---|---
Probability of | 85.043% | 85.042% | 85.04% | 85.02% | 84.95% | 84.47% | 82.63% | 79.27%
interaction (%) |  |  |  |  |  |  |  |
Ratio to | 1 | 0.99998 | 0.99993 | 0.99973 | 0.99893 | 0.99320 | 0.97163 | 0.93215
on-centre beam |  |  |  |  |  |  |  |
Table 2: Probability of interaction of protons hitting the target at different
offset from the target center

###  III.3 Code development

For coding convenience, the target is divided into two volumes, upstream and
downstream, with the break at MC0 of horn 1 (see Fig. 5). In this initial
study, the code only allows target position to change in integer steps of
sphere diameter. To locate the target at any arbitrary position, the code
needs to be able to split a sphere into two volume at MC0, using volume
subtraction method.

Figure 5: HepRApp visualizaion of horn 1 (inner conductor) and spherical array
target

In LBNE g4lbne source code, new methods for upstream and downstream physical
volume placement, namely
LBNEVolumeplacements::PlaceFinalUpstrMultiSphereTarget and
LBNEVolumePlacements::PlaceFinalDownstrMultiSphereTarget, were written in
parallel with the volume placement methods for simple cylinder and simple box.

In the original code, the number of graphite fin targets and target position
with respect to horn 1 is calculated based on surveyed information entered by
user. This will allow users to run the simulation based on actual
measurements. LBNEVolumePlacements::SegmentTarget() was branched out do a
similar task on SAT. The method is able to adjust target length and position
to get an integer number of spheres in total, in upstream and downstream
sections.

The above methods are called after an user interface (UI) card
fUseMultiSphereTarget. Another new UI cards in LBNEPlacementMessenger was also
defined to take sphere radius as an input. These cards are used in macros
files to change the configurations.

The new geometry was verified using plots of Geantino tracks, HepRApp and
OpenGL visualization. Once the new geometry has been verified to work,
simulations were run on Fermilab computer clusters. Each configuration was run
for 500 jobs, each with 100,000 Proton on Target (POT) at beam power 1.2MW.
Jobs of the same configuration were merged together and make comparison
histograms of neutrino flux at the Far Detector. The neutrino of interest are
unoscillated muon neutrino $\nu_{\mu}$, electron neutrino $\nu_{e}$ and their
anti-partners compared to LBNE Nominal design.

Figure 6: OpenGL visualization of target and inner conductor of horn 1 with a
proton on target resulting in an interaction.

##  IV Result and Analysis

The target was pushed inside horn 1 deeper, using less than 3mm radial
clearance (which is the clearance used by LBNE Nominal design). Using target
diameter of 17mm and two interaction length of material (969mm) Fig. 7 shows a
higher neutrino flux at 2-3 GeV except for the most retracted target (shown in
gray). The higher flux is possibly due to the geometry of the SAT with more
than two interaction length at center.

The ratio plot in Fig. 8 shows the gain ratio of flux compared to LBNE Nominal
graphite target. The violet plot has a similar target length inside horn 1. To
measure the first and second oscillation maximum, the range of interest is 0.5
- 4 GeV. A more favorable flux is shown by the red plot, where the target is
pushed deeper inside and flux gain is suppressed at high energy. A sharp
transition as the target is retracted from the deepest position (shown in red)
by only 51mm (blue plot) begs for more in depth investigation.

Figure 7: Neutrino flux produced with SAT at the Far Detector. Sphere
diameter: 17mm, beam sigma: 17/6mm, target length: 969mm. Legend indicates
length of target outside and inside horn 1 in mm, respectively. Nominal plot
(black) refers to flux produced by LBNE Graphite target with length inside
horn 1 is approximately 504 mm, comparable to the violet plot. Figure 8: Ratio
of neutrino flux produced with SAT over LBNE Nominal target. Sphere diameter:
17mm, beam sigma: 17/6mm, target length: 969mm.

Plot in Fig. 9 shows a relatively stable neutrino flux for the chosen beam
size of R/3. Flux loss is approximately 1% when beam offset is 1mm. This
result is roughly in agreement with prediction of probability of proton
interaction in target in Table. 2.

Figure 9: Ratio of neutrino flux produced by SAT at the Far Detector of off-
centre beam over on-centre beam. Sphere diameter: 17mm, beam sigma: 17/6mm,
target length: 969mm. Figure 10: Ratio of neutrino flux produced with 13mm SAT
at the Far Detector of off-centre beam over LBNE Nominal target. Sphere
diameter: 13mm, beam sigma: 13/6mm, target length: 975mm

Decreasing sphere diameter to 13mm (which elongates the target to 975mm)
shifts the flux gain to higher energy and introduce flux loss at lower energy,
shown in Fig. 10. In most target positions, flux is loss at the first
oscillation maximum of 0.8GeV. Flux increases up to 50% at 6 GeV, which is not
of much interest at 1300km.

Figure 11: Ratio of neutrino flux produced with SAT at the Far Detector of
off-centre beam over LBNE Nominal target. Sphere diameter: 17mm, beam sigma:
1.7mm, target length: 969mm

A smaller beam size of 1.7 mm was also used on 17mm diameter spheres for
simulation, result is shown in Fig. 11. The flux gain is higher at low energy,
compare to the beam size of 2.8mm. A possible explanation is when the beam is
more concentrated in the center with longer path length, produced particles
traverse through more material and as a consequence, loose energy and produce
lower energy neutrino. The plot also shows that, the less target extension
outside horn 1, the better flux gain at 0-3.5GeV. This is better demonstrated
by the plot in Fig. 12. With the same target length inside horn 1 (shown in
red), a shorter target give a better flux scheme: higher gain at 0-3.5GeV, and
more flux suppression at after 3.5 GeV.

Figure 12: Ratio of neutrino flux produced with SAT with less than two
interaction lengths at the Far Detector of off-centre beam over LBNE Nominal
target. Sphere diameter: 17mm, beam sigma: 1.7mm, target length: 901mm

##  V conclusion

The spherical array target with less than two interaction length and a smaller
beam sigma of 1.7mm gives higher flux at 0-3.5 GeV and suppresses flux at
energy higher than 3.5 GeV. For spherical array target, the beam is relatively
stable with less than 2% flux loss as the beam is offset from target center by
1mm.

A further study of spherical array target should isolate the effect of
different parameters and study each separately: material (graphite vs.
Beryllium), geometry (fins vs. spherical), total target length and target
length inside/outside horn 1.

More engineering study is required to obtain accurate size of Helium can, size
and material of helical spacer around target spheres and stress of spherical
array target under a beam size of 1.7mm.

The source code also need to be upgraded to move the target to any arbitrarily
position and to include the helical spacer. It’s also possible to show number
of protons interacted in target and interaction location to better understand
the effect of target size and geometry.

##  VI acknowledgement

I am grateful to my advisors, Dr. Alberto Marchionni , Dr. Paul Lebrun and Dr.
Laura Fields at Fermilab for their support and guidance.

## References

  * (1) Milind Diwan, Robert Svoboda, and James Strait. The long-baseline neutrino experiment. European Strategy Preparatory Group 2012.
  * (2) LBNE Collaboration. Long-Baseline Neutrino Experiment (LBNE) Project Conceptual Design Report Volume 2: The Beamline at the Near Site. Fermi National Accelerator Laboratory 2012.
  * (3) CJ Densham, O Caretta, TR Davenne, MD Fitton, P Loveridge, M Rooney. Conceptual Design Study of the LBNE Target and Beam Window, Final Report. STFC Rutherford Appleton Laboratory  2010
  * (4) T. Nakadaira. Performance and operational feedback of T2K graphite target. T2K Collaboration 2012