\section{Socialist Science}
\label{sec:SocialistScience}

Drug company Q reads a scholarly article in a distinguished journal authored by academic A.  Motivated in significant part by A's findings, Q starts a new research effort in this seemingly promising direction.

One year and one million dollars later, Q realizes A's findings are basically wrong~\cite{begley2012drug,prinz2011believe,freedman2015economics,chalmers2009avoidable,scott2008design,begley2013unappreciated,steckler2015preclinical,de2015failure,kola2004can,macleod2014biomedical,kyzas2007almost,hirst2014need,Hyman155cm11,miller2010pharma,schoenfeld2013everything}.

Variants of this story may be invoked by a convenient buzz phrase:  science's ``replication crisis,'' ``replicability crisis,'' or ``reproducibility crisis.''  Q is not always a drug company.  A is usually an academic.

Q, the sucker in this story, realizes he is continually getting shafted~\cite{ioannidis2005most,ioannidis2013s,pankevich2014improving,lowenstein2009uncertainty,henderson2013threats,dirnagl2012international,dirnagl2009stroke,dirnagl2006bench,rosenblatt2016incentive}.  Systematic attempts to repeat many published studies fail~\cite{mobley2013survey,open2015estimating,ioannidis2005contradicted,ioannidis2007non,steward2012replication,chang2015economics}, spreading awareness to the masses and stoking discontent.  Debates erupt over the scale of the problem~\cite{maxwell2015psychology,anderson2016response,stroebe2014alleged,klein2014investigating,pashler2012replicability,camerer2016evaluating,etz2016bayesian,michalek2010costs}.  Some hope extermination of psychology departments will be sufficient~\cite{fanelli2010positive,lilienfeld2012public,cesario2014priming,Pashler01112012,bakker2012rules,Bones01052012,wagenmakers2011psychologists,ferguson2015everybody}.  Others work to re(de)fine the notion of ``replication''~\cite{goodman2016does,brandt2014replication,clemens2015meaning}.  Incremental improvements to the current system are proposed~\cite{ioannidis2014make,collins2014nih,vesterinen2010improving,asendorpf2013recommendations,nosek2015promoting,cumming2013new,iqbal2016reproducible,holdren2013increasing,miguel2014promoting,roche2014troubleshooting,chalmers2014increase,ioannidis2014increasing,salman2014increasing,yordanov2015avoidable,moher2015four,ioannidis2014assessing,vasilevsky2013reproducibility,international2004international,kilkenny2010improving,macleod2009reprint,tooth2005quality,festing2002guidelines,moher2010guidance,casadevall2012reforming,nosek2012scientificI,nosek2012scientificII,everett2015tragedy,nekrutenko2012next,sandve2013ten,ioannidis2011improving,pusztai2013reproducibility,valentine2011replication,kidwell2016badges,koole2012rewarding,landis2012call,moher2010consort,moseley2014beyond,peers2012search,peers2014can,ram2013git,schooler2014metascience}~\footnote{{\em{Open science}}, or {\em{even-more-socialist science}}, is a common theme among proposed incremental improvements.  The idea of encouraging more transparent access to data and analysis code is an attractive one.  We ourselves pushed it in particle physics -- hard, and for many years.  Unfortunately, if unsurprisingly, incentives are simply too misaligned for it to work.  If we, as a society, want more of something -- like apples, say, or knowledge about how nature works -- we may be better off making it easy for people who produce apples to sell them than mandating that they make their orchards available for anyone to pick.}.  Most are never enacted~\cite{baker2014two}.  The rest have limited impact~\cite{vanpaemel2015we,roche2015public,vines2014availability,wicherts2006poor,moher2016increasing,grant2013reporting,joseph2013open,song2010dissemination,prayle2012compliance,clarke2010clinical,florez2016bias,bramhall2015quality,plint2006does}.

The incentive structures contributing to the replication crisis are observed to be deeply entrenched features of the ecosystem within which science is done~\cite{begley2015reproducibility,simmons2011false,fiedler2011voodoo,ferguson2012vast,young2008current,tsilidis2013evaluation,brembs2013deep,steen2013has,oboyle2014chrysalis,stamatakis2013undue,sena2010publication,mathieu2009comparison,hannink2013comparison,crowe2015patients,blumle2016fate,chan2014increasing,glasziou2014reducing,glasziou2014role,hoffmann2012scatter,franco2014publication,button2013power,dwan2008systematic,ioannidis2008most,van2010can,ioannidis2012science,stroebe2012scientific,wagenmakers2012agenda,john2012measuring,fanelli2011negative,alberts2014rescuing,ioannidis2014publication,chambers2014instead} and subsequently applied~\cite{glasziou2005paths,duff2010adequacy,dancey2010quality,kilkenny2009survey,mcglynn2003quality,lemon2005surveying,ioannidis2001completeness,savovic2012influence,turner2008selective,bero1998closing,ramagopalan2014prevalence,williamson2012driving}.  This ecosystem comprises institutions (including funding agencies, industry, universities, and academic journals), culture, accepted practice (including procedures by which grant money is awarded, articles published, promotions granted, and accolades bestowed), and the bureaucracy required to support all of these.  The ecosystem is complex, multifaceted, and not easily changed.  The incentive structures contributing to the replication crisis are, similarly, not easily changed.  A timely solution to the replication crisis probably requires a new ecosystem.

For convenience, we refer to the current ecosystem~\cite{stephan1996economics} -- funded by taxpayers, with published results available to all citizens -- as {\em{socialist science}}~\footnote{The phrase ``socialist science'' is obviously a grotesquely crude caricature of the intricate and often nuanced set of incentives joining the actors in the current science ecosystem.  We intend the phrase as a neutral description of an aspect of the current ecosystem germane to the present discussion.  A reader who dislikes the phrase is encouraged to mentally replace it with ``the ecosystem within which science is currently carried out,'' or some alternative reference thereto.}~\footnote{We focus on the incentive flaws common to all socialist science, ignoring differences among scientific disciplines.  Rather than treat the symptoms, which express differently in the social sciences, life sciences, and physical sciences, we focus on the underlying disease.}.

\section{Capitalist Science}
\label{sec:CapitalistScience}

Perhaps we can construct a different ecosystem.

Lacking imagination, we look to see what other professions do.  Many professions exhibit a very interesting behavior.  There occur {\em{pairwise transactions}}, in which one party {\em{sells}} something of value to another party.  One party (A) gives something of value to some other party (Q).  In return, Q gives A something called {\em{money}}.

Perhaps our academic (A) can sell what he learns to those (Q) who find A's information valuable.  Many deep pocketed industries, including the pharmaceutical industry, make high stakes decisions motivated by academic research.  A should have little difficulty finding customers Q.

To incentivize accuracy, we need some sort of transaction that causes A to lose money if he turns out to be wrong.  We therefore need transactions for which the (in)accuracy of A's information can eventually be objectively determined, and we need a reliable, low-cost procedure for making that determination.

Suppose Q has a specific question and is willing to pay for a useful answer.  Q can ensure the usefulness of any answer he receives by specifying, explicitly or categorically, the set of allowed possible answers.  A third party (X) brokers the transaction, holds money from Q and A in escrow, and acts as arbiter.  If A provides an answer outside Q's set of allowed possible answers, A loses money.

In many cases of practical interest, Q can be tasked with determining the accuracy of A's answer.  If Q wants to know what chemical matter binds to a particular target and A answers with a specific molecule, Q will eventually know whether A's answer was right or wrong, and Q can provide evidence to this effect to arbiter X.  X collects a deposit from Q at the beginning of the transaction, and X returns this deposit to Q when Q provides sufficient evidence by a previously agreed upon date.  Q pays the same amount whether A's answer turns out to be right or wrong.  X's fee for brokering the transaction is independent of whether A's answer turns out to be right or wrong and whether or not Q provides evidence X considers sufficient.  A has an explicit monetary incentive to be accurate and an explicit monetary disincentive to provide information unless he is pretty sure he is correct.  Q has an explicit monetary incentive to ask a question that will be directly useful to him, and for which he can eventually determine, with reasonable proof, the accuracy of whatever answer A provides.  X, as an ongoing business interest, has a clear monetary incentive to thoughtfully serve its role as a reasonable and objective adjudicator of the evidence Q provides.  Q, A, and X all have an interest in ensuring Q's question is carefully specified to avoid subsequent ambiguities as to A's correctness.  This transaction protocol enforces a set of incentives facilitating paid, arm's length transfer of useful, bluntly honest information from A to Q.

Using this transaction protocol, introduced for the first time in Ref.~\cite{knuteson2016knowledge} and implemented at Ref.~\footnote{Kn-X; \url{http://kn-x.com}.  Patent pending.}, scientists can sell what they learn from their research.

For convenience, we refer to this new ecosystem as {\em{capitalist science}}~\footnote{The phrase ``capitalist science'' is intended as a neutral description of a salient feature of this new ecosystem.  A reader who dislikes the phrase is encouraged to mentally replace it with ``the new ecosystem,'' or something similar.}~\footnote{The ``capitalist science'' in this article supersedes that of Ref.~\cite{knuteson2011capitalist}, which in retrospect is more of a hybrid between socialist and capitalist science.}~\footnote{To be clear, our use of the term ``capitalist'' does not arise from a belief, seemingly held by many, that free markets are the optimal solution to all problems, nor from a belief that global financial markets at the time of this writing function well and should be emulated.  (You have no idea.)  The complexity and frequent opacity of today's capitalism highlights the glaring need for a mechanism facilitating useful, bluntly honest information transfer between remote parties.   Given the embarrassing, hidden-in-plain-sight, farcically tragic comedy of errors that is recent financial history, a mechanism facilitating useful, bluntly honest, arm's length information transfer may turn out to be our best shot at saving capitalism \ldots\ or at least at sending a few of those responsible to prison next time around.}.

In capitalist science, A makes money only if he turns out to be correct.  A loses money if he turns out to be wrong.  Every transaction includes a monetary incentive for Q to determine the accuracy of A's answer and to back this determination with evidence deemed sufficient by an objective third party (X).  These features directly address socialist science's reproducibility crisis.  These incentives, present in capitalist science, explicitly and directly reward accuracy and penalize inaccuracy.  The absence of such explicit and direct incentives in socialist science is the reason socialist science is suffering a reproducibility crisis.

Capitalist science is the solution to socialist science's reproducibility crisis.

\section{Summary}
\label{sec:Summary}

Now that the solution to socialist science's reproducibility crisis has been found, we can roll up our sleeves \ldots\ sit back, relax, and let events unfold.

Socialist science will continue much as it has, producing results of similar quality.  Serious, thoughtful, and impassioned attempts will be made to save it.  Some of these may appear to work for a while, but ultimately they will fail.  The reproducibility crisis has been stewing for decades~\cite{lykken1968statistical,elms1975crisis,greenwald1975consequences,rosenthal1979file,altman1994scandal,hackam2006translation,pocock1987statistical,sterling1995publication,vul2009puzzlingly,easterbrook1991publication,kerr1998harking}.  The relevant malincentives are too firmly embedded in its large and unwieldy ecosystem for socialist science to produce a meaningful solution from within.

The key to change is Q, the socialist science sucker from Section~\ref{sec:SocialistScience}.  Q is smart.  Q is a socialist science sucker only because socialist science fails to sufficiently discourage the publication of inaccurate information by failing to hold A sufficiently accountable for being wrong.  Q is a socialist science sucker only because Q has had no other option.

Q now has another option.  Before launching a one year, one million dollar research project motivated by a result produced by socialist science, Q can spend ten days and ten thousand dollars on a few questions to check it out.  Before investing ten million dollars on a new materials science innovation, venture capitalists can anonymously conduct due diligence with a convenience and at a depth hitherto unimaginable.  Relying on free information of questionable accuracy produced by socialist science can be very costly.  Q is smart; Q is accustomed to paying for things of value; and Q's stakes are high.

Over time, Q increasingly relies on information obtained from capitalist science, where A has skin in the game.  When A has truly valuable information, he finds himself more interested in selling it than publishing it.  For less robust results, A is more inclined to publish, avoiding the harsh penalty imposed by capitalist science for providing information that turns out to be wrong.  Over time, useful, robust information slowly migrates to capitalist science, where it is appropriately valued.  Socialist science continues to publish the results A is willing to give away for free.

The forces in capitalist science's favor are so strong that socialist science should be allowed to fade into irrelevance in its natural course~\footnote{Although the information market unleashed by capitalist science could create millions of new science-related jobs, it would be irresponsible to reduce funding to socialist science until that promise has been realized.}.  Socialist science has been a remarkable institution.  The knowledge it has provided has been extraordinary in power and in scope.  Its use of an incentive structure spurned by the most developed of today's modern economies only makes its achievements all the more impressive.

We fully hope and expect socialist science will linger among us for many years.  At some deep level, the purpose and goals of socialist science are fundamentally different from those of capitalist science.  Capitalist science is designed to facilitate useful, accurate information transfer, and to provide A with a strong incentive to learn information some Q will find valuable.  Socialist science, in contrast, is not.


\section{Introduction}
It has been well known that disorder may change the behavior of the system,
not only quantitatively but also qualitatively, in stochastic systems,
as well as in classical or quantum spin systems and the quantum localization
problem\cite{RBS,IM}. The effect of disorder is more remarkable when the
spatial dimension of the system is lower. For example, at one-particle level,
a localization where a particle stays at a certain specific site with much
higher probability than at the other sites\cite{BG}, and a phase separation,
where the regions with high particle density and low particle density are
separated spontaneously\cite{Krug}, are well known phenomena.

In many cases, the disorder which is considered is a random disorder:
A disorder is represented as identically and independently distributed
(i.i.d.) stochastic variables and thus, by definition,
these variables are not spatially
correlated each other. In a real situation, however, disorder may be
spatially correlated and it is valuable to consider a model with spatially
correlated disorder.

Aperiodic disorder is a class of disorder with spatial correlation.
An aperiodic disorder is generated deterministically  from a set of specific
rules. That is the point which distinguishes an aperiodic disorder from
random disorder and which enables us to study systems with aperiodic
disorder systematically. Therefore study of systems with aperiodic disorder
is a good first step toward understanding of more general system with
spatially correlated disorder. Although some investigations have been carried
out about the effect of aperiodic disorder in classical and quantum spin
chains\cite{Luck,Hermisson,Vieira,Hida} so far, to our knowledge,
just a few in stochastic
systems. It has been known that even at one-particle level, an interesting
behavior can be observed: An anomalous diffusion where the variance grows
less slowly than linear growth in time \cite{ITR} and correspondingly the
stationary probability distribution with multifractality\cite{Miki}.
Note that systems with aperiodic disorder are not only of theoretical and
mathematical interest, but also they have been artificially fabricated and
investigated experimentally\cite{apexp}.

One of the next interesting problems is to make clear an interplay between
quenched disorder with spatial correlation and interaction between particles.
For this purpose we study a simple multiparticle stochastic process,
called the zero-range process (ZRP)\cite{EH} with aperiodic disorder.
In ZRP, it is remarkable that the exact stationary measure can be written
in a factorized form, which is derived from the one-particle stationary
probability distribution. Moreover, under a certain condition, condensation
can be observed, where a great part of the particles condense at one
specific site. The effect of random quenched disorder in ZRP has already
been investigated\cite{JSI,GL}. In this paper, specifically, we consider
ZRP with disorder constructed from an aperiodic sequence,
called the paperfolding (PF) sequence\cite{Dekking}.
This is one of the models where a multifractal stationary probability
distribution can be observed.

\section{A zero-range process with aperiodic disorder}
Let us consider a one-dimensional lattice with periodic boundary condition.
Each site $j=1,2,\cdots, L$ can accommodate any number of particles.
A particle hops from site to site with designated hop rates and only hopping
to one of the nearest neighboring sites is allowed. The hop rates are quenched
variables, where we denote the forward rate from site $j$ to $j+1$ by $p_j$
and the backward rate from $j$ to $j-1$ by $q_{j-1}$ (see Fig.\ref{zrp}).
\begin{figure}
\begin{center}
\includegraphics[width=8cm]{ZRP.eps}
\caption{ZRP on one-dimensional lattice with aperiodic disorder.
\label{zrp}}
\end{center}
\end{figure}
Here we assume that the hop rates do not depend on the number of particles
at the departure site (sometimes called the "pure chipping process"\cite{LMS}),
although generally the hop rates in ZRP are given by a function of the number
of particles at the departure site. For example, in a noninteracting particle
system the hop rates are directly proportional to the number of particles
at the departure site.

Let us construct the aperiodic disorder according to the PF sequence.
The PF sequence $S=AABAABBA\cdots$ is generated systematically
by the initial sequence $S_1=AA$ and the substitution rules,
$AA \rightarrow AABA$, $AB \rightarrow AABB$, $BA \rightarrow ABBA$, and
$BB \rightarrow ABBB$. The length of the sequence of the $n$-th generation,
$S_n$, is $2^n$ and in the $n \rightarrow \infty$ limit the ratio of the number
of A's, $N_n(A)$, to that of B's, $N_n(B)$, converges to unity.
Then the hop rates are assigned as
\begin{equation}
p_j=1,\quad \text{for all $j$},
\label{forwardrates}
\end{equation}
and
\begin{equation}
q_j=
\begin{cases}
a, &\quad \text{the $j$-th symbol of $S$ is $A$},
\\
b, &\quad \text{the $j$-th symbol of $S$ is $B$}.
\end{cases}
\label{backwardrates}
\end{equation}

At one-particle level, the drift velocity $v_d$ is proportional to the
difference between unity and the product of $q_j/p_j$ \cite{Derrida}
\begin{equation}
v_d \propto  1-\prod_{j=1}^L \frac{q_j}{p_j}.
\label{driftvel}
\end{equation}

The state of the multiparticle system is specified by the particle
configuration
$\{n_j\}=\{n_1,n_1,\cdots,n_L\}$, where $n_j$ denotes the number of particles
at site $j$. The total number of particles, $N=\sum_{j=1}^Ln_j$ is conserved,
since the periodic boundary condition is imposed.

Let $P(\{n_j\};t)$ denote the probability that the configuration $\{n_j\}$
is observed at time $t$. The master equation is given as
\begin{eqnarray}
\frac{\partial}{\partial t} P(\{n_j\};t)
&=& \sum_j p_{j-1} P(\{\cdots, n_{j-1}+1,n_{j}-1,\cdots \};t)
\nonumber
\\
&&+\sum_j q_{j} P(\{\cdots, n_{j}-1,n_{j+1}+1,\cdots \};t)
\nonumber
\\
&&-\sum_j (p_{j}+q_{j-1})P(\{n_j\};t).
\label{mastereq}
\end{eqnarray}
One of the remarkable features of ZRP is that the exact multiparticle
stationary measure  $P(\{n_j\})$ can be obtained,
from the one-particle distribution $\{P_j\}$
(the stationary probability that the particle exist at the
$j$-th site), in a factorized form:
\begin{equation}
P(\{n_j\})=\frac{1}{Z_{L,N}}\prod_{j=1}^L(P_j)^{n_j}\delta(\sum_{j=1}^Ln_j,N),
\label{prodmeasure}
\end{equation}
where $Z_{L,N}$ denotes the partition function
\begin{equation}
Z_{L,N}=\sum_{\{n_j\}}\prod_{j=1}^L(P_j)^{n_j}\delta(\sum_{j=1}^Ln_j,N).
\label{partitionfunction}
\end{equation}
The symbol $\delta(\sum_{j=1}^Ln_j,N)$ indicates that the total number of
particles is fixed at $N$.

The one-particle distribution $\{P_j\}$ is given by the stationary solution
of the master equation describing the random walk on the lattice:
\begin{equation}
\frac{\partial P_j(t)}{\partial t}
= p_{j-1}P_{j-1}(t) + q_{j}P_{j+1}(t)
-(p_{j}+q_{j-1})P_j(t).
\label{randomwalk}
\end{equation}
It is known that the stationary solution can be obtained
analytically\cite{Derrida}:
\begin{equation}
P_j = \frac{r_j}{\sum_{k=1}^Lr_k},
\label{probdist}
\end{equation}
where
\begin{equation}
r_j = \frac{1}{p_j}\left(1+
\sum_{i=1}^{L-1}\prod_{k=1}^i\frac{q_{j+k-1}}{p_{j+k}}\right).
\label{rj_aux}
\end{equation}

Under the condition for the drift velocity Eq.(\ref{driftvel}) to vanish
\begin{eqnarray}
b&=&a^{-\delta(n)},
\nonumber
\\
&&\delta(n)=N_n(A)/N_n(B),
\label{zerocurrent}
\end{eqnarray}
the stationary distribution $\{P_j\}$ shows multifractality and therefore
the distribution is neither extended nor localized\cite{Miki}.
Hereafter we restrict ourselves to this case. This model has only one
free parameter $a$.

\section{Condensation and scaling of the number of particles out of condensate}
Figure \ref{profile} shows the density profile obtained from Monte Carlo
simulations.
\begin{figure}
\begin{minipage}{0.8\hsize}
\begin{center}
\includegraphics[width=8cm]{lowdensity.eps}
\end{center}
\end{minipage}
\begin{minipage}{0.8\hsize}
\begin{center}
\includegraphics[width=8cm]{highdensity.eps}
\end{center}
\end{minipage}
\caption{Density profile obtained by simulations. (a) Low-density case,
$L=256$ and $N=16$. (b) High-density case, $L=N=256$. The inset is the
same plot with the ordinate in the log-scale.
\label{profile}}
\end{figure}
When the average density is low, the interaction between particles is
ineffective and the profile shows a hierarchical structure, as observed for
the one-particle distribution (see Fig.\ref{profile}(a)).
On the other hand, when the average density increases, many particles
condense at the site which gives the maximum of the one-particle
distribution (see Fig.\ref{profile}(b)).
This condensation is already known to occur for randomly disordered case.
An analogy to the Bose-Einstein condensation is pointed out\cite{EH}.

In order to study the condensation in more detail, we introduce the
grand-canonical(GC) formulation. In the GC formulation, the grand partition
function ${\mathcal Z}_L(z)$ is given in a factorized form:
\begin{eqnarray}
{\mathcal Z}_L(z) &=& \sum_{N \ge 0} z^N Z_{L,N}
\nonumber
\\
&=& \prod_{j=1}^L \frac{1}{1-zP_j},
\label{grandpartition}
\end{eqnarray}
where $z$ denotes the fugacity. Then the local density at site $j$ is
obtained as
\begin{equation}
\rho^{\rm GC}_j(z) = \frac{zP_j}{1-zP_j}.
\label{localdensity}
\end{equation}
The average density is given as
\begin{equation}
\rho^{\rm GC}(z) = \frac{1}{L}\sum_{j=1}^L\frac{zP_j}{1-zP_j}.
\label{avdensity}
\end{equation}
From these expressions, it is immediately understood that the density of
the site, at which the maximum of the one-particle distribution is given,
rapidly increases, since at the site, as $z$ increases, the numerator of
Eq.(\ref{localdensity}) increases most rapidly and the denominator
decreases most rapidly. However, for a finite system,
Eqs.(\ref{localdensity}) and (\ref{avdensity}) can take any finite value
for $0 \le z <1/P_{\max}$. Therefore, the GC formulation cannot be broken
as long as we consider a finite system. As is well known, condensation
is characterized by the breaking of the GC formulation. Then it is necessary
to consider the thermodynamic limit where $L,N \rightarrow \infty$ with
$N/L$ fixed to be finite.

In the thermodynamic limit, the expression of the total density
Eq.(\ref{avdensity}) is rewritten as
\begin{eqnarray}
\rho^{\rm GC}(z) &=& \int_0^1dx\frac{zP(x)}{1-zP(x)},
\qquad(x:=j/L, L \rightarrow \infty)
\nonumber
\\
&=& \int_0^{1/\max\{P_j\}}dpf(p)\frac{zp}{1-zp},
\label{TDdensity}
\end{eqnarray}
where the summation in the r.h.s. of Eq.(\ref{avdensity}) is replaced with
the integral, and $f(p)$ denotes the distribution function of $p$.
A condensation can occur if the integral converges in the limit
$z \rightarrow 1/\max\{P_j\}$. For the integral to converge, the distribution
function $f(p)$ must be an increasing function of $(p_{\max}-p)$ in the
neighborhood of $p_{\max}=\max\{P_j\}$. In the case where the one-particle
distribution $\{P_j\}$ has multifractality, this condition is considered
to be satisfied. However, it is still quite
difficult to calculate the integral and evaluate the critical density at
which the condensation occurs. This consideration is similar to that
for the density of states near the lowest energy level in Bose-Einstein
condensation. Hereafter, we consider the case where the total density
is large enough for the condensation to occur in the thermodynamic
limit.

The total number of particles out of the condensate, $N_{\rm out}$, is
evaluated as
\begin{eqnarray}
N_{\rm out} &=& \sum_{j \ne j_{\max}} \rho_j^{\rm GC}(z)
\nonumber
\\
&=& \sum_{j \ne j_{\max}} \frac{zP_j}{1-zP_j}
\nonumber
\\
&\approx& \sum_{j \ne j_{\max}} \frac{P_j}{P_{\max}-P_j},
\label{nouteq}
\end{eqnarray}
where $j_{\max}$ is the condensation site and in the last line the fugacity
$z$ is replaced with $1/P_{\max}$.

Figure \ref{nout} shows the dependence of $N_{\rm out}$
on the system size $L$ for $a=0.3$, $0.4$, and $0.5$. It is found that
$N_{\rm out}$ increases algebraically with respect to $L$:
\begin{equation}
N_{\rm out} \sim L^{\gamma(a)},
\label{noutpower}
\end{equation}
where the power-law exponent $\gamma$ depends on $a$. This result is
qualitatively different from that in the case with random disorder,
where $N_{\rm out}$ is kept ${\mathcal O}(1)$, independent of the
system size\cite{JSI}. This can be considered as one of the
characteristic results generated by the spatially
correlated aperiodic disorder and interaction between particles.
\begin{figure}
\begin{center}
\includegraphics[width=8cm]{noutplot.eps}
\caption{System size dependence of the number of particles out of the
condensate, $N_{\rm out}(L)$, for $a=0.3$, $0.4$, and $0.5$. The lines are
guide for the eyes.
\label{nout}}
\end{center}
\end{figure}

Figure \ref{gma} shows the dependence of the exponent
$\gamma$ on $a$. It is a monotonically increasing function: The stronger the
disorder, the smaller the exponent. However, we are not sure about the values
to which the exponent converges in the $a \rightarrow 0$ (extremely strong
disorder) and $a \rightarrow 1$ (extremely weak disorder) limits. It is quite
difficult to obtain them numerically, since for $a \rightarrow 0$ some
measures are quite small, and for $a \rightarrow 1$ the difference of the
measures is quite small.

\begin{figure}
\begin{center}
\includegraphics[width=8cm]{gamma_a.eps}
\caption{Dependence of power-law exponent on disorder strength,
$\gamma(a)$ in Eq.(\ref{noutpower}) for the PF model and the BBP.
\label{gma}}
\end{center}
\end{figure}

In our previous study, we found that the multifractal structure of
the one-particle distribution can be well reproduced by the binomial
branching process (BBP)\cite{Miki}. This process is constructed by
the iteration of dividing the segments into two halves and
uneven partitioning of the measure by assiging one segment to $p$ and the
other to $(1-p)$\cite{CMKRS}, where $p$, assumed to be $>1/2$,
is the only free parameter of the process.
The disorder strength of the PF model, $a$, and the parameter $p$ can be
related through the multifractal singular exponent:
\begin{equation}
\alpha_{\min}(a) = -\log_2p.
\label{pf_bbp}
\end{equation}
For the distribution generated by the BBP with $p$, which corresponds to
the given disorder strength $a$, we evaluate $N_{\rm out}$ through
Eq.(\ref{nouteq}). The evaluated $N_{\rm out}$ also shows a power-law
dependence on the system size $L$,
similar to that for the PF model. Moreover, the
exponent $\gamma$ agrees very well, as shown in Fig.\ref{gma}.
The one-particle distributions themselves for these two processes are
not so similar each other,
although their multifractal structures are similar. This suggests that
the multifractal structure of the one-particle distributionis essential.
It should be noted that for the PF model the structure is resulted from
the aperiodic disorder.

\section{Summary}
We investigated a phenomena emerged by both an aperiodic disorder and
interaction between particles using the ZRP with disorder according to
the PF sequence, for which the underlying one-particle stationary
distribution is multifractal. In the stationary state, when the density
is sufficiently high, a condensation occurs, where most particle
condense at one specific site. For a finite size system, the number of
particles out of the condensate increases with the system size
algebraically, contrary to the randomly disordered case where it is
suppressed to ${\mathcal O}(1)$. The distribution generated by the BBP
has a similar multifractality and can reproduce the exponent well.
This is a characteristic result by an interplay between the effects of
spatially correlated disorder and that of the interaction between
particles.

\begin{acknowledgments}
This research was supported by the initiative-based project E-05
"Creation and Sustainable Governance of New Commons through Foundation of
Integrated Local Environmental Knowledge (ILEK)", Research Institute for
Humanity and Nature (RIHN).
\end{acknowledgments}


\section{Introduction}
\label{sec:intro}
Let
\begin{equation}\label{eq:sampling_densities}
 X\sim N_p(\theta,\sigma^2I_p),\ S\sim \sigma^2\chi^2_n,
\end{equation}
where $X$ and $S$ are independent and $\theta$ and $\sigma^2$ are both unknown,
and where
\begin{equation}\label{eq:dim}
 p\geq 3, \quad n\geq 3.
\end{equation}
Consider the problem of estimating the mean vector $\theta$
under the loss function
\begin{equation}\label{eq:scaled_quadratic_loss}
 L(\{\theta,\sigma^2\},d)=\|d-\theta\|^2/\sigma^2.
\end{equation}
We study the question of admissibility/inadmissibility of shrinkage-type estimators
of the form
\begin{equation}\label{eq:delta_phi}
 \delta_\phi(X,S)=\left(1-\phi(W)/W\right)X,
\end{equation}
where $ W=\|X\|^2/S$.
We do so by examining the existence of solutions to a differential inequality
which arises from an unbiased estimate of the difference in risk
between $\delta_\phi$ and $\delta_{\phi+g}$.
Hence we are more properly studying what may be termed quasi-admissibility
and quasi-inadmissibility of such estimators.
Quasi-inadmissibility implies inadmissibility under conditions of risk finiteness,
while quasi-admissibility is relatively weaker.

Stein in his unpublished lecture notes, \cite{Brown-1988}, \cite{Bock-1988}, \cite{Rukhin-1995} and
\cite{Brown-Zhao-2009} among others have studied the admissibility
question from this point of view (without necessarily using the term quasi-admissibility)
under known $\sigma^2$.
Of course, \cite{Brown-1971} has largely settled the admissibility/inadmissibility question
when $\sigma^2$ is known.

Our efforts focus generally on finding a boundary between quasi-admissibility and
quasi-inadmissibility for shrinkage estimators of the form \eqref{eq:delta_phi}.
(See Theorem \ref{thm:main}) We also apply the result to a class of generalized Bayes
estimators related to the class of \cite{Strawderman-1971} priors for the known variance
problem and establish a boundary on the tail behavior which also separates
quasi-admissibility from quasi-inadmissibility.

While minimaxity of shrinkage estimators in the unknown scale case has been extensively
studied by many authors, relatively little is known about admissibility in this case.
\cite{Strawderman-1973} and \cite{Zinodiny-etal-2011}
gave a class of proper Bayes minimax and hence
admissible estimators under unknown $\sigma^2$.
Note that proper Bayes estimators by \cite{Strawderman-1973} and \cite{Zinodiny-etal-2011}
are not of the form given by \eqref{eq:delta_phi} whereas
generalized Bayes estimators by \cite{Maruyama-2003a},
\cite{Maru-Straw-2005} and \cite{Maruyama-Strawderman-2009} are of this form.

While our results on quasi-admissibility do not settle the admissibility issue,
it seems likely to us that generalized Bayes estimators satisfying our conditions
for quasi-admissibility are admissible, perhaps under mild additional conditions.
We are decidedly not claiming that such a result would be easily established!
Certainly those found to be quasi-inadmissible are in fact inadmissible under conditions
of finiteness of risk.

An unbiased estimator of of the risk, $R(\{\theta,\sigma^2\},\delta_\phi)$,
for an estimator of the form
\eqref{eq:delta_phi} is given by
\begin{equation}\label{eq:unbiased_estimator_1}
 p+(n+2)D_\phi(W)
\end{equation}
where
\begin{equation}\label{eq:unbiased_estimator_D_phi}
D_\phi(w)=\frac{\{\phi(w)-2c_{p,n}\}\phi(w)}{w}-d_n \phi'(w)\left\{1+\phi(w)\right\},
\end{equation}
with
\begin{equation}\label{eq:unbiased_estimator_c_d}
 c_{p,n}=(p-2)/(n+2)\text{ and }d_n=4/(n+2).
\end{equation}
This result follows from
\citeapos{Stein-1981} identity and well known identities for chi-square distributions
(see e.g.~\cite{Efron-Morris-1976}).
We may refer to \eqref{eq:unbiased_estimator_1} as a SURE estimate of risk and
to \eqref{eq:n+2D} below as a SURE estimate of difference in risk.
A sufficient condition for its validity is
that $\phi$ be absolutely continuous
and that each term of $E[D_\phi(W)]$ be finite.
Let $\varPhi$ be a family of functions $\phi$, satisfying these sufficient conditions,
\begin{equation}\label{conditions_varphi}
 \varPhi =\left\{\phi: E[D_\phi(W)]<\infty, \text{ absolute continuous}\right\}.
\end{equation}

If $\delta_{\phi+g}$ is of the form \eqref{eq:delta_phi} with $\phi(w)$ replaced by
$\phi(w)+g(w)$, an unbiased estimator of the difference in risk between $\delta_\phi$
and $\delta_{\phi+g}$,
\begin{equation*}
 R(\{\theta,\sigma^2\},\delta_\phi)-R(\{\theta,\sigma^2\},\delta_{\phi+g})
\end{equation*}
is given by
\begin{equation}\label{eq:n+2D}
\begin{split}
 (n+2)\Delta(w;\phi,g)&=(n+2)\{D_\phi(w)-D_{\phi+g}(w)\} \\
&=(n+2)g(w)\{\Delta_1(w;\phi)+\Delta_2(w;\phi,g)\}
\end{split}
\end{equation}
where
\begin{align}\label{Delta_1}
\Delta_1(w;\phi)=2\frac{c_{p,n}-\phi(w)}{w}+d_n\phi'(w)
\end{align}
and
\begin{align}\label{Delta_2}
\Delta_2(w;\phi,g)=\frac{-g(w)}{w}+d_ng'(w)+d_n \frac{g'(w)}{g(w)}\{1+\phi(w)\}.
\end{align}
One may find an estimator dominating $\delta_\phi$
by finding a non-zero solution $g(\cdot)\in\varPhi$ to the
differential inequality $\Delta(w;\phi,g)\geq 0$, where $\Delta(w;\phi,g)$ is given by \eqref{eq:n+2D},
providing the resulting estimator has finite risk.
Here is the definition of quasi-admissibility and quasi-inadmissibility used in this paper:
\begin{definition}\label{def:quasi}
 \begin{enumerate}
  \item An estimator $\delta_\phi$ of the form \eqref{eq:delta_phi} is said to be
quasi-admissible if any solution $g(w)\in \varPhi$ of the inequality
$\Delta(w;\phi,g)\geq 0$ satisfies $g(w)\equiv 0$,
\item $\delta_\phi$ is said to be
quasi-inadmissible if there exists a solution, $g(w)\in \varPhi$, which
is non-vanishing on some open interval, to the differential inequality $\Delta(w;\phi,g)\geq 0$.
 \end{enumerate}
\end{definition}
For technical reasons we will restrict the class of $\phi(\cdot)$ to the subclass
$\varPhi_{A}$ of $\varPhi$, defined as follows,
\begin{align}
 \varPhi_{A}=\left\{\phi\in\Phi, \text{and }\phi \text{ satisfies \ref{AA1}, \ref{AA2}, \ref{AA3}, and \ref{AA4} below}\right\},
\end{align}
\begin{enumerate}[label= A\arabic*]
\item \label{AA1} $\phi(0)=0$ and $\phi(w)\geq 0$ for any $w\geq 0$,
\item \label{AA2} $\phi(w)$ has at most finitely many local extrema,
\item \label{AA3} $\phi'(w)$ has only finitely many discontinuities and $\phi'(w)$ is continuous
from the right at $0$.
\item \label{AA4} $\liminf_{w\to\infty} w\phi'(w)/\phi(w)\geq 0$ and
       $\limsup_{w\to\infty} w\phi'(w)/\phi(w)\leq 1$.
\end{enumerate}

Note that James-Stein-type estimators $(1-a/W)X$ with $\phi(w)\equiv a$ do not satisfy Assumption \ref{AA1}.
However such estimators are inadmissible
and are dominated by the positive part version $(1-a/W)_+ X$ for which $\phi_+(w)=\min(w,a)$.
$\phi_+(w)$ does in fact satisfy Assumption \ref{AA1}. The positive part modification of any $\delta_\phi$
for which $\lim_{w\to 0}\phi(w)>0$ will similarly satisfy $\phi_+(0)=0$.
Assumption \ref{AA2} assumes that $\phi(w)$ does not oscillate excessively and that
$\lim_{w\to\infty}\phi(w)$ exists. Assumption \ref{AA3} is
used in controlling the local behavior of $\phi$ and of $\phi'$.
Assumptions \ref{AA1}--\ref{AA4} are
satisfied by linear estimators of the form $\delta(X)=\alpha X$
for $0\leq \alpha \leq 1$ and for which $\phi(w)=(1-\alpha) w$.
These estimators are unique proper Bayes and admissible in the normal case for $0\leq\alpha<1$.
As far as we know, Assumptions \ref{AA1}--\ref{AA4} cover the positive part version
of all minimax estimators in the literature.
We emphasize that while we address quasi-admissibility and inadmissibility only for
$\delta_\phi$ for $\phi\in\Phi_A\subset\Phi$, we allow competitive estimators
of the form $\delta_{\phi + g}$ for $g\in\Phi$.

In Section \ref{sec:main} we will show the following result, which establishes
\begin{equation*}
  \phi(w)= \frac{p-2}{n+2}- \frac{\beta_\star}{\log w}
\end{equation*}
as the asymptotic boundary between quasi-admissibility and quasi-inadmissibility
where
\begin{equation}\label{beta_star_intro}
 \beta_\star=\frac{d_n(1+c_{p,n})}{2}=\frac{2(p+n)}{(n+2)^2}.
\end{equation}
\begin{description}
	 \item[Quasi-admissibility:]
If $\phi\in\varPhi_A$ and there exists $w_*$ and $b<1$
such that
\begin{align}\label{eq:DD1.intro}
 \phi(w)\geq \frac{p-2}{n+2}- b\frac{\beta_\star}{\log w}, \ \forall w\geq w_*,
\end{align}
		    then $\delta_\phi$ is quasi-admissible.
\item[Quasi-inadmissibility:]
If $\phi\in\varPhi_A$ and there exists $w_*$ and $b>1$
such that
\begin{align}\label{eq:DD2.intro}
 \phi(w)\leq \frac{p-2}{n+2}-b\frac{\beta_\star}{\log w}, \ \forall w\geq w_*,
\end{align}
then $\delta_\phi$ is quasi-inadmissible (and hence inadmissible).
\end{description}
In Section \ref{sec:Bayes}, we find a generalized Bayes estimator with
asymptotic behavior
\begin{equation*}
\lim_{w\to\infty} \log w\left(\frac{p-2}{n+2}-\phi(w)\right)=b\beta_\star,
\end{equation*}
for all $b>0$. The corresponding generalized prior is given by
\begin{equation*}
\frac{1}{\sigma^2}\times\frac{1}{\sigma^p}G(\|\theta\|/\sigma)
\end{equation*}
with
\begin{equation*}
 G(\|\mu\|)=
\int_0^1 \left(\frac{\lambda}{1-\lambda}\right)^{p/2}
\exp\left(-\frac{\lambda}{1-\lambda}\frac{\|\mu\|^2}{2}\right)\lambda^{-2}
\left(\log\frac{1}{\lambda}\right)^b
  d\lambda.
\end{equation*}
Hence, $b<1$ and $b>1$ imply quasi-admissibility and quasi-inadmissibility, respectively,
of the associated generalized Bayes estimators.
Interestingly,
the boundary $b=1$ also appears in the known $\sigma^2$ case when estimating $\mu$
with $Z\sim N_p(\mu,I_p)$. By using \citeapos{Brown-1971}
sufficient condition, the generalized Bayes estimator with respect to $G(\|\mu\|)$ above
is admissible (resp.~inadmissible) when $b\leq 1$ (resp.~$b>1$).
This nice correspondence leads naturally to the conjecture:
a quasi-admissible generalized Bayes estimator satisfying \eqref{eq:DD1.intro} is admissible.

An extension to the general class of spherically symmetric distributions is briefly
considered in Section \ref{sec:ssd}. We give some concluding remarks in Section \ref{sec:cr}.
Some technical proofs are given in \ref{sec:AP}ppendix.

 \section{Quasi-admissibility}
\label{sec:main}
The main result of this paper, Theorem \ref{thm:main}, gives sufficient conditions for quasi-admissibility
and quasi-inadmissibility for estimators $\delta_\phi$ of the form \eqref{eq:delta_phi},
for $\phi\in\varPhi_A$.
In preparation, we first give several lemmas.
Recall that the unbiased estimator of the difference in risk between $\delta_\phi$
and $\delta_{\phi+g}$ is given by
\begin{equation}\label{eq:n+2D_sec2}
\begin{split}
 (n+2)\Delta(w;\phi,g)&=(n+2)\{D_\phi(w)-D_{\phi+g}(w)\} \\
&=(n+2)g(w)\{\Delta_1(w;\phi)+\Delta_2(w;\phi,g)\}
\end{split}
\end{equation}
where
\begin{align}\label{Delta_1_sec2}
\Delta_1(w;\phi)=2\frac{c_{p,n}-\phi(w)}{w}+d_n\phi'(w)
\end{align}
and
\begin{align}\label{Delta_2_sec2}
\Delta_2(w;\phi,g)=\frac{-g(w)}{w}+d_ng'(w)+d_n \frac{g'(w)}{g(w)}\{1+\phi(w)\},
\end{align}
and where $c_{p,n}=(p-2)/(n+2)$ and $d_n=4/(n+2)$.
Note that $\Delta_2(w;\phi,g)$ is well-defined for $w$ such that $g(w)\neq 0$,
but $\Delta(w;\phi,g)$ is well-defined even when $g(w)=0$.

The first lemma gives necessary conditions on $g(w)$ for $\Delta(w;\phi,g)$
to be nonnegative for all $w\geq 0$.
\begin{lemma}\label{lem:GG}
Suppose $\Delta(w;\phi,g)\geq 0$ for all $w\geq 0$ with $\phi\in\varPhi_A$ and $g\in\varPhi$.
Then
\begin{enumerate}
\item \label{BB1} $g(0)\geq 0$,
\item \label{BB2} $ g(w)\geq 0 $ for all $w> 0$,
\item \label{BB3} Suppose $g(w_0)>0$. Then, for any $w\geq w_0$, $g(w)>0$.
\end{enumerate}
\end{lemma}
\begin{proof}
 Section \ref{sec.ap.1} in Appendix.
\end{proof}
Recall that finiteness of $E[\phi(W)^2/W]$ is a necessary condition for $\phi$ to be in $\varPhi$.
Lemma \ref{CC0} below provides a necessary condition for
$E[\phi(W)^2/W]$ to be finite and hence for $\phi$ to be in $\varPhi$.
It is needed in the proof of Lemma \ref{lem:CC}.
\begin{lemma}\label{CC0}
A necessary condition for $E[\phi(W)^2/W]$ to be finite for any $(\theta,\sigma^2)$ is that
\begin{equation}\label{eq:liminfphi}
 \liminf_{t\to\infty}|\phi(t)|^{d_n}/t=0.
\end{equation}
\end{lemma}
\begin{proof}
 Section \ref{sec.ap.2} in Appendix.
\end{proof}
 Let $\mathcal{G}\subset \varPhi$ be a class of nonnegative functions which satisfy \ref{BB1}, \ref{BB2} and
 \ref{BB3} of Lemma \ref{lem:GG} and Lemma \ref{CC0}.
 The following lemma is key in proving the main result.
Recall that Assumption \ref{AA2} assumes that $\phi(w)$ does not oscillate excessively and that
$\lim_{w\to\infty}\phi(w)$ exists.
In the following lemma, let $\phi_*=\lim_{w\to\infty} \phi(w) \in [0,\infty]$ and
\begin{equation}\label{beta_star_star}
 \beta_\star=\frac{d_n(1+c_{p,n})}{2}=\frac{2(p+n)}{(n+2)^2}.
\end{equation}
\begin{lemma}\label{lem:CC}
 Suppose $\phi\in\varPhi_A$.
 \begin{enumerate
\item \label{CC1}
      Suppose $\phi_*<\infty$ and there exists $w_0$ and $b<1$ such that
\begin{align}\label{eq:CC1.lem}
 \phi(w)\geq \frac{p-2}{n+2}- b\frac{\beta_\star}{\log w}, \ \forall w\geq w_0.
\end{align}
\begin{enumerate}[label= (1.\alph*), ref= 1.\alph*]
 \item \label{CC1.1}For all $w\geq w_0$,
\begin{equation}\label{eq:CC1.1.lem}
 \Delta_1(w;\phi)-d_n\phi'(w)-\frac{2b\beta_\star}{w\log w}\leq 0.
\end{equation}
\item \label{CC1.2}For any $g\in\mathcal{G}$ except $g\equiv 0$,
there exists $w_*\in(w_0,\infty)$ such that
\begin{equation}\label{eq:CC1.2.lem}
\Delta_2(w_*;\phi,g)+d_n\phi'(w_*)+\frac{2b\beta_\star}{w_*\log w_*}<0.
\end{equation}
\end{enumerate}
\item \label{CC2}
      Suppose $\phi_*=\infty$.
\begin{enumerate}[label= (2.\alph*), ref= 2.\alph*]
 \item \label{CC2.1}
       There exists $w_0$ such that
\begin{equation}\label{eq:CC2.1.lem}
 \Delta_1(w;\phi)+d_n\phi'(w)\leq 0,\text{ for all }w\geq w_0.
\end{equation}
 \item  \label{CC2.2}
	For any $g\in\mathcal{G}$ except $g\equiv 0$,
there exists $w_*\in(w_0,\infty)$ such that
\begin{equation}\label{eq:CC2.2.lem}
\Delta_2(w_*;\phi,g)-d_n\phi'(w_*)<0.
\end{equation}
\end{enumerate}
 \item \label{CC3}
Suppose there exists $w_0$ and $b>1$ such that
\begin{align}\label{eq:CC3.lem}
 \phi(w)\leq \frac{p-2}{n+2}-b\frac{\beta_\star}{\log w}, \ \forall w\geq w_0.
\end{align}
\begin{enumerate}[label= (3.\alph*), ref= 3.\alph*]
 \item \label{CC3.1}
       There exists $w_1$ such that
       \begin{equation}\label{eq:CC3.1.lem}
 \Delta_1(w;\phi)-\frac{2b\beta_\star}{w\log w}\geq 0,\text{ for all }w\geq w_1.
       \end{equation}
 \item \label{CC3.2}
       Fix
\begin{equation}
 \nu=\min\left(1,\frac{2b\beta_\star-d_n(1+\phi_*)}{2d_n(3+\phi_*)}\right).
\end{equation}
       Let $k(w)\in \mathcal{G}$ be any non-decreasing continuous function with
$k(0)=0$, $w_\sharp=\sup\{w: k(w)=0\}$  and $k(\infty)=1$.
Then there exists $w_*$, independent of $k(w)$, such that
\begin{equation}\label{eq:CC3.2.lem}
 \Delta_2(w;\phi,k(w)\{\log (w+e)\}^{-1-\nu})

+\frac{2b\beta_\star}{w\log w} \geq 0
\end{equation}
       for all $w> \max(\max(w_*,w_1), w_\sharp)$ and  $e=\exp(1)$.
 \end{enumerate}
\end{enumerate}
\end{lemma}
\begin{proof
Section \ref{sec.ap.3} in Appendix.
\end{proof}
Note, in part \ref{CC3}, $ \Delta_2(w;\phi,k(w)\{\log (w+e)\}^{-1-\nu})$ is well-defined
for $w>w_\sharp$ by the definition of $\Delta_2$ given by \eqref{Delta_2_sec2}.

\smallskip

The following result is the main result of this section.
\begin{thm}\label{thm:main}
 Suppose $\phi\in\varPhi_A$.
 \begin{enumerate
  \item \label{DD1}
[quasi-admissibility] If there exists $w_*$ and $b<1$
such that
\begin{align}\label{eq:DD1}
 \phi(w)\geq \frac{p-2}{n+2}- b\frac{\beta_\star}{\log w}, \ \forall w\geq w_*,
\end{align}
then $\delta_\phi$ is quasi-admissible.
  \item \label{DD2}
	[quasi-inadmissibility]
If there exists $w_*$ and $b>1$
such that
\begin{align}\label{eq:DD2}
 \phi(w)\leq \frac{p-2}{n+2}-b\frac{\beta_\star}{\log w}, \ \forall w\geq w_*,
\end{align}
then $\delta_\phi$ is quasi-inadmissible (and hence inadmissible).
 \end{enumerate}
\end{thm}
\begin{proof}\mbox{}
 [Part \ref{DD1}]
 By Parts \ref{CC1} ($\phi_*<\infty$) and \ref{CC2} ($\phi_*=\infty$) of Lemma \ref{lem:CC},
there exists $w_*$ such that $ \Delta_1(w_\star;\phi)+\Delta_2(w_\star;\phi,g)<0$
for any $g\in\mathcal{G}$ except $g\equiv 0$.
 Therefore any solution $g(w)\in\mathcal{G}$ of the differential inequality
\begin{equation*}
 g(w)\left\{\Delta_1(w;\phi)+\Delta_2(w;\phi,g)\right\}\geq 0
\end{equation*}
must be identically equal to $0$,
or equivalently $\delta_\phi$ is quasi-admissible.

 [Part \ref{DD2}]
By \eqref{eq:DD2}, we have $\phi_*\leq (p-2)/(n+2)=c_{p,n}$ and hence
 \begin{equation}\label{futoushiki_2.1}
  d_n(1+\phi_*)\leq 2\beta_\star<2b\beta_\star
 \end{equation}
since $b>1$.
As in Part \ref{CC3} of Lemma \ref{lem:CC}, let
\begin{equation}
 \nu=\min\left(1,\frac{2b\beta_\star-d_n(1+\phi_*)}{2d_n(3+\phi_*)}\right).
\end{equation}
 Take any $k(w)$ with $ w_\sharp=\max(w_1, w_*)$ where $w_1$ and $w_1$ are both determined by
Part \ref{CC3} of Lemma \ref{lem:CC}.
Let  $g(w)=k(w)\{\log (w+e)\}^{-1-\nu}\in\mathcal{G}$.
 Then we have
\begin{align*}
 \Delta(w)= g(w)\left\{\Delta_1(w;\phi)+\Delta_2(w;\phi,g)\right\}
 \begin{cases}
 =0 & 0\leq w\leq w_\sharp \\
\geq 0 & w> w_\sharp,
 \end{cases}
\end{align*}
 where $ \Delta(w)=0$ for $0\leq w\leq w_\sharp$ since $g(w)=0$ and $ \Delta(w)\geq 0$ for $w> w_\sharp$
 since
\begin{equation*}
 \left( \Delta_1(w;\phi)-\frac{2b\beta_\star}{w\log w} \right)
  +\left(\Delta_2(w;\phi,g)

+\frac{2b\beta_\star}{w\log w}\right)\geq 0
\end{equation*}
by Part \ref{CC3} of Lemma \ref{lem:CC}.
Hence $\delta_\phi$ is quasi-inadmissible.
\end{proof}

\begin{remark}
 Note that it is possible that an estimator which is quasi-admissible
according to the above definition may fail to be admissible for several reasons.
Here are some of them.
First, there may be an estimator that is not of the form  \eqref{eq:delta_phi}
that dominates $\delta_\phi$. Second,
there may be an estimator of the form  \eqref{eq:delta_phi} with $g(w)\notin \varPhi$
that dominates $\delta_\phi$.
Third there may be an estimator that dominates $\delta_\phi$
but does not satisfy the differential inequality $\Delta(w;\phi,g)\geq 0$.
Hence quasi-admissibility is quite weak as an optimality criterion.

Quasi-inadmissibility, on the other hand, is more compelling in the sense that
if $\delta_\phi$ is quasi-inadmissible then it is inadmissible and dominated by $\delta_{\phi+g}$.
Note that requiring both $\phi$ and $g$ to be in $\varPhi$ implies that
the risk of $\delta_{\phi+g}$ is finite.
\end{remark}

\subsection{General spherically Symmetric distributions}
\label{sec:ssd}
We may also study the more general canonical spherically symmetric setting where
$(X,U)$ has a spherically symmetric density of the form
\begin{equation}\label{density_f.ssd}
 \sigma^{-p-n}f(\{\|x-\theta\|^2+\|u\|^2\}/\sigma^2).
\end{equation}
Here the $p$-dimensional vector $X$ has mean vector $\theta$, the $n$-dimensional ``residual''
vector $U$ has mean vector $0$ and $(X,S)$ is sufficient, where $S=\|U\|^2$.
The scale parameter, $\sigma^2$, is assumed unknown.
Consider the problem of estimating the mean vector $\theta$
under the loss function
\begin{equation}\label{eq:scaled_quadratic_loss.ssd}
 L(\{\theta,\sigma^2\},d)=\|d-\theta\|^2/\sigma^2.
\end{equation}
The most important such setting is the Gaussian case
\begin{equation}\label{eq:sampling_densities.ssd}
 X\sim N_p(\theta,\sigma^2I_p),\ S\sim \sigma^2\chi^2_n,
\end{equation}
which is studied in Section \ref{sec:main}, but there is considerable interest
in the case of heavier tailed distributions such as the multivariate-$t$.

In the general spherically symmetric case,
\eqref{eq:unbiased_estimator_1} is not an unbiased estimate of risk but has been used as a substitute for such an estimator.
In particular, if $(X,S)$ has density \eqref{density_f.ssd} and $F(\cdot)$ is defined as
\begin{equation}
 F(t)=\frac{1}{2}\int_t^\infty f(v)dv.
\end{equation}
Then as essentially shown by several authors in various settings
(see e.g.~\cite{Kubokawa-Srivastava-2001} and \cite{Fourdrinier-Strawderman-2014})
\begin{align*}
&R(\{\theta,\sigma^2\},\delta_{\phi}) \\ &=p+(n+2)\int_{\mathbb{R}^{p+n}}D_\phi(w)
 \frac{F(\{\|x-\theta\|^2+\|u\|^2\}/\sigma^2)}{\sigma^{p+n}}dx du,
\end{align*}
where $ D_\phi(w)$ is given in \eqref{eq:unbiased_estimator_D_phi}.
Hence, in this setting,
 \begin{align*}
&  R(\{\theta,\sigma^2\},\delta_\phi)-R(\{\theta,\sigma^2\},\delta_{\phi+g}) \\
  &=(n+2)\int_{\mathbb{R}^{p+n}} g(w)\{\Delta_1(w;\phi)+\Delta_2(w;\phi,g)\}
  \frac{F(\{\|x-\theta\|^2+\|u\|^2\}/\sigma^2)}{\sigma^{p+n}}dx du
 \end{align*}
 where $w=\|x\|^2/\|u\|^2$.
Thus, study of existence of solutions to $\Delta(w)\geq 0$ is relevant in the general spherically
symmetric case as well as in the Gaussian case,
and defining quasi-admissibility/inadmissibility as in Definition \ref{def:quasi} implies
that Theorem \ref{thm:main} remains valid in this more general setting.

\section{Generalized Bayes estimators in the Normal case}
\label{sec:Bayes}
\subsection{Known variance case}
Let $Z\sim N_p(\mu,I_p)$. Consider estimation of $\mu$ under quadratic loss $\|\hat{\mu}-\mu\|^2$.
The MLE, $Z$ itself, is inadmissible for $p\geq 3$ as shown in \cite{Stein-1956}.
\cite{Brown-1971} showed that admissible estimators should be proper Bayes or generalized Bayes
estimators with respect to an improper prior and gave
a sufficient condition for generalized Bayes estimators to be admissible or inadmissible.

Let the prior be of the form
\begin{equation}\label{general_spherical_prior_known}
 \pi(\mu)=G(\|\mu\|;a,L)
\end{equation}
where
\begin{equation}\label{GGG}
 G(\|\mu\|;a,L)=
\int_0^1 \left\{\frac{\lambda}{1-\lambda}\right\}^{p/2}
  \exp\left(-\frac{\lambda}{1-\lambda}\frac{\|\mu\|^2}{2}\right)
  \lambda^a L(1/\lambda)  d\lambda
\end{equation}
where $p/2+a+1>0$. We assume the following on $L$:$[1,\infty)$ $\to$ $[0,\infty)$
\begin{enumerate}[label= L\arabic*]
 \item \label{L.A.3}
       $L(y)$ is slowly varying at infinity, that is, for all $c>0$,
\begin{equation*}
 \lim_{y\to\infty}L(cy)/L(y)=1.
\end{equation*}
 \item \label{L.A.1}
       $L(y)$ is ultimately monotone,
 \item \label{L.A.2}
       $L(y)$ is differentiable with ultimately monotone derivative $ L'(y)$,
\end{enumerate}
By Proposition 1.7 (11) of \cite{geluk-dehaan-1987}, Assumptions \ref{L.A.3} and \ref{L.A.2} implies
\begin{equation*}
 \lim_{y\to \infty}y\frac{L'(y)}{L(y)}=0.
\end{equation*}
Under the prior given by \eqref{GGG}, the marginal density is
\begin{align*}
 m(\|z\|;a,L)&=\int_{\mathbb{R}^p}\frac{1}{(2\pi)^{p/2}}
 \exp\left(-\frac{\|z-\mu\|^2}{2}\right)G(\|\mu\|;a,L)d\mu \\
&=\int_0^1
 \exp\left(-\frac{\lambda\|z\|^2}{2}\right)\lambda^{p/2+a}L(1/\lambda)d\lambda \\
&=\int_0^\infty
 \exp\left(-\frac{\lambda\|z\|^2}{2}\right)f(\lambda;a,L)d\lambda,
\end{align*}
where $f(\lambda;a,L)=\lambda^{p/2+a}L(1/\lambda)I_{(0,1)}(\lambda)$.
Note that $f(\lambda;a,L)$ is ultimately monotone as a function of $1/\lambda$ since
\begin{enumerate}
 \item When $p/2+a=0$, $L(1/\lambda)$ itself is ultimately monotone.
 \item When $p/2+a\neq 0$, $ \lim_{\lambda\to 0}\lambda f'(\lambda)/f(\lambda)=p/2+a\neq 0$.
\end{enumerate}
Since $f(\lambda;a,L)$ is ultimately monotone and since
$m(\|z\|;a,L)$ is the Laplace transform of $f$,
a Tauberian Theorem
(See e.g.~\cite{Feller-1971} Theorem 13.5.4) implies that
\begin{equation}\label{TaubTaub}
 m(\|z\|;a,L)  \approx \Gamma(p/2+a+1)\left(2/\|z\|^{2}\right)^{p/2+a+1}L(\|z\|^2)
\end{equation}
as $\|z\|\to\infty$.
As shown in Appendix \ref{sec:ap.bound},
$\|z\|\|\nabla \log m(\|z\|;a,L)\|$ is bounded.
By Theorem 6.4.2 of \cite{Brown-1971},
divergence (convergence) of the integral
\begin{equation}\label{BrownBrown}
 \int_1^\infty \frac{dr}{r^{p-1}m(r;a,L)}
\end{equation}
corresponds to admissibility (inadmissibility) of a generalized Bayes estimator
with bounded $\|z\|\|\nabla \log m(\|z\|;a,L)\|$.
Hence, by \eqref{TaubTaub} and \eqref{BrownBrown},
divergence (convergence) of the integral
\begin{equation}\label{BrownBrown}
 \int_1^\infty \frac{r^{2a+3}}{L(r^2)}dr
\end{equation}
corresponds to admissibility (inadmissibility).
It is clear that $a>-2$ and $a<-2$ imply admissibility and inadmissibility, respectively.

When $a=-2$, the fact that
\begin{equation}
 \int_1^\infty \frac{dr}{r\{\log r\}^b}
  \begin{cases}
   = \infty & b\leq 1 \\ <\infty & b>1
  \end{cases}
\end{equation}
is helpful to determine the boundary.
Since
\begin{equation}\label{eq:log_b}
\int_1^\infty \frac{r^{2a+3}}{L(r^2)}dr
=
\frac{1}{2^b}\int_1^\infty \frac{1}{r\{\log r\}^b}\frac{\{\log r^2\}^b}{L(r^2)}dr,
\end{equation}
we have a following result on admissibility and inadmissibility of the
(generalized) Bayes estimator with respect to $G(\|\mu\|;a,L)$.
\begin{thm}[Admissibility]\label{thm:ad}\mbox{}
\begin{enumerate}
 \item \label{thm:ad.p1}
       Suppose $a>-2$. The generalized Bayes estimator is inadmissible.
 \item \label{thm:ad.p2}
       Suppose $a=-2$ and $\log(y)/L(y)$ is ultimately monotone non-decreasing.
       The generalized Bayes estimator is admissible.
\end{enumerate}
\end{thm}
\begin{thm}[Indmissibility]\label{thm:inad}\mbox{}
\begin{enumerate}
 \item \label{thm:inad.p1}Suppose $a=-2$ and $\{\log(y)\}^b/L(y)$ for $b>1$ is ultimately monotone non-increasing.
The generalized Bayes estimator is inadmissible.
 \item \label{thm:inad.p2}
       Suppose $a<-2$. The generalized Bayes estimator is inadmissible.
\end{enumerate}
\end{thm}
\begin{remark}[A boundary esimator for the known variance case]\label{rem:b.1}\mbox{}
 Consider the particular choice
\begin{equation}\label{eq:Lb}
a=-2\text{ and } L(1/\lambda)=\left(\log\frac{1}{\lambda}\right)^b \text{ for }b>0.
\end{equation}
Then the prior is given by
\begin{equation*}
\int_0^1 \left\{\frac{\lambda}{1-\lambda}\right\}^{p/2}
  \exp\left(-\frac{\lambda}{1-\lambda}\frac{\|\mu\|^2}{2}\right)
  \lambda^{-2} \left(\log\frac{1}{\lambda}\right)^b  d\lambda.
\end{equation*}
By following \cite{Strawderman-1971}, the corresponding generalized Bayes estimator
is $(1-\psi_{-2,b}(\|Z\|^2)/\|Z\|^2)Z$ where
\begin{equation*}
 \psi_{-2,b}(v)=v\frac
  {\int_0^1\lambda^{p/2-1}\{\log(1/\lambda)\}^b\exp(-v\lambda/2)d\lambda}
  {\int_0^1\lambda^{p/2-2}\{\log(1/\lambda)\}^b\exp(-v\lambda/2)d\lambda},
\end{equation*}
As shown in Appendix \ref{sec.ap.5}, we have
\begin{equation}\label{known_variance_phi_asymptotic}
 \lim_{v\to\infty}\left(\log v\right)\left\{p-2-\psi_{-2,b}(v)\right\}=2b.
\end{equation}
 Hence
 by Part \ref{thm:ad.p2} of Theorem \ref{thm:ad} and
Part \ref{thm:inad.p1} of Theorem \ref{thm:inad},
the generalized Bayes estimator with asymptotic behavior
 \begin{equation*}
  \left(1-\left\{p-2 - \frac{b}{\log \|Z\|}\right\}\frac{1}{\|Z\|^2}\right)Z
 \end{equation*}
is admissible and inadmissible for $ b\leq 1$ and $ b>1$.
Thus the estimator
 \begin{equation*}
  \left(1-\left\{p-2 - \frac{1}{\log \|Z\|}\right\}\frac{1}{\|Z\|^2}\right)Z
 \end{equation*}
 is a boundary estimator. See also Corollary 6.3.2 of \cite{Brown-1971} and
 Theorem 6.1.1 of \cite{Strawderman-Cohen-1971} for related discussions, but where
 the $b/\log \|z\|$ term is not included.
\end{remark}

\subsection{Unknown variance case}
Let $X$ and $S$ be given by \eqref{eq:sampling_densities} and let the prior be of the form
\begin{equation}\label{general_spherical_prior}
 \pi(\theta,\sigma^2)=\frac{1}{\sigma^2}\pi(\theta|\sigma^2)=\frac{1}{\sigma^2} \times \frac{1}{\sigma^p}G(\|\theta\|/\sigma)
\end{equation}
where $G$ is given by \eqref{GGG} and $1/\sigma^2$ is a standard non-informative prior for $\sigma^2$.

The following two theorems relate quasi-admissibility/inadmissibility in the unknown
variance case to admissibility/inadmissibility in the known variance case as given in
Theorems \ref{thm:ad} and \ref{thm:inad}.
\begin{thm}[Quasi-admissibility]\label{thm:q.ad}\mbox{}
 \begin{enumerate}
 \item Suppose $a>-2$. The generalized Bayes estimator is quasi-admissible.
 \item Suppose $a=-2$ and $\{\log(y)\}^b/L(y)$ for $b<1$ is monotone non-decreasing.
       The generalized Bayes estimator is quasi-admissible.
 \end{enumerate}
\end{thm}
\begin{thm}[Quasi-indmissibility]\label{thm:q.inad}\mbox{}
\begin{enumerate}
 \item Suppose $a=-2$ and $\{\log(y)\}^b/L(y)$ for $b>1$ is monotone non-increasing.
       The generalized Bayes estimator is quasi-inadmissible.
 \item Suppose $a<-2$. The generalized Bayes estimator is quasi-inadmissible.
\end{enumerate}
\end{thm}
\begin{proof}[Proof of Theorems \ref{thm:q.ad} and \ref{thm:q.inad}]
By following \cite{Maru-Straw-2005} and \cite{Maruyama-Strawderman-2009},
the generalized Bayes estimator under the prior given by \eqref{general_spherical_prior}
is $\delta_\phi$ with
\begin{equation*}
 \phi_{a,L}(w)
  =w\frac{\int_0^1\lambda^{p/2+a+1}L(1/\lambda)(1+w\lambda)^{-(p+n)/2-1}d\lambda}
  {\int_0^1\lambda^{p/2+a}L(1/\lambda)(1+w\lambda)^{-(p+n)/2-1}d\lambda}.
\end{equation*}
By a change of variables ($t=w\lambda$), we have
\begin{equation*}
  \phi_{a,L}(w)
=\frac{\int_0^w t^{p/2+a+1}L(w/t)(1+t)^{-(p+n)/2-1}dt}
  {\int_0^w\lambda^{p/2+a}L(w/t)(1+t)^{-(p+n)/2-1}dt}.
\end{equation*}
By Assumption \ref{L.A.3} and the Lebesgue dominated convergence theorem,
\begin{equation}\label{phi.L.limit}
 \begin{split}
 \lim_w \phi_{a,L}(w)
&=\frac{\int_0^\infty t^{p/2+a+1}(1+t)^{-(p+n)/2-1}dt}
  {\int_0^\infty t^{p/2+a}(1+t)^{-(p+n)/2-1}dt} \\
&=\frac{p/2+a+1}{n/2-a-1}
\end{split}
\end{equation}
which is increasing in $a$ and is equal to $(p-2)/(n+2)$ when $a=-2$.
Hence, by Theorem \ref{thm:main},
$a>-2$ and $a<-2$ implies quasi-admissibility and quasi-inadmissibility, respectively.

When $a=-2$, take
\begin{equation}
 L(1/\lambda)=\left\{\log\frac{1}{\lambda}\right\}^{b}
\end{equation}
for $b>0$ and consider
\begin{equation}\label{eq:phi-2b}
 \phi_{-2,b}(w)
  =w\frac{\int_0^1\lambda^{p/2-1}\{\log(1/\lambda)\}^b(1+w\lambda)^{-(p+n)/2-1}d\lambda}
  {\int_0^1\lambda^{p/2-2}\{\log(1/\lambda)\}^b(1+w\lambda)^{-(p+n)/2-1}d\lambda}.
\end{equation}
Then we have
\begin{equation}\label{lim.final.qa}
 \lim_{w\to\infty}(\log w)\left(\frac{p-2}{n+2}- \phi_{-2,b}(w)\right)=b\frac{2(p+n)}{(n+2)^2}=b\beta_\star
\end{equation}
where $\beta_\star$ is given by Theorem \ref{thm:main}.
See Section \ref{sec.ap.5} in Appendix for the derivation of \eqref{lim.final.qa}.
Further
the inequality
\begin{equation}\label{final.ineq}
 \begin{split}
 \phi_{-2,b}(w)&=\frac{\int_0^1\lambda^{p/2-1}\{\log(1/\lambda)\}^b(1+w\lambda)^{-(p+n)/2-1}d\lambda}
  {\int_0^1\lambda^{p/2-2}\{\log(1/\lambda)\}^b(1+w\lambda)^{-(p+n)/2-1}d\lambda} \\
 &=\frac{\displaystyle \int_0^1\lambda \frac{\{\log(1/\lambda)\}^b}{L(1/\lambda)}
 \lambda^{p/2-2}L(1/\lambda)(1+w\lambda)^{-(p+n)/2-1}d\lambda}
 {\displaystyle \int_0^1 \frac{\{\log(1/\lambda)\}^b}{L(1/\lambda)}
 \lambda^{p/2-2}L(1/\lambda)(1+w\lambda)^{-(p+n)/2-1}d\lambda} \\
 & \leq (\geq)
 \frac{ \int_0^1\lambda
 \lambda^{p/2-2}L(1/\lambda)(1+w\lambda)^{-(p+n)/2-1}d\lambda}
 {\int_0^1
 \lambda^{p/2-2}L(1/\lambda)(1+w\lambda)^{-(p+n)/2-1}d\lambda} \\ &=\phi_{-2,L}(w)
\end{split}
\end{equation}
follows for $b>0$
when $ \{\log(y)\}^b/L(y)$ is monotone non-increasing (non-decreasing).
From \eqref{phi.L.limit}, \eqref{lim.final.qa}, \eqref{final.ineq} and Theorem \ref{thm:main},
 the two theorems follow.
\end{proof}
\begin{remark}[A boundary esimator for the unknown variance case]\label{rem:b.2}\mbox{}
 For the unknown variance case, Theorem \ref{thm:main} established the boundary
 between quasi-admissibility and quasi-inadmissibility for estimator
 of the form $(1-\phi(W)/W)X$ as
 \begin{equation}\label{eq:bound.unknown}
  \left(1-\left\{\frac{p-2}{n+2} - \frac{b\beta_*}{\log \|W\|}\right\}\frac{1}{\|W\|^2}\right)X
 \end{equation}
 with $b<1$ corresponding to quasi-admissibility and
 $b>1$ corresponding to quasi-inadmissibility.
The generalized prior
\begin{equation*}
 \pi(\theta,\sigma^2)=\frac{1}{\sigma^2} \times \frac{1}{\sigma^p}G(\|\theta\|/\sigma)
\end{equation*}
with $ G$ given by \eqref{GGG} where
\begin{equation}\label{eq:Lb.1}
a=-2\text{ and } L(1/\lambda)=\left(\log\frac{1}{\lambda}\right)^b \text{ for }b>0
\end{equation}
leads to a generalized Bayes estimator with $\phi$ given in \eqref{eq:phi-2b}.
As shown in Appendix, the asymptotic behavior of this $\phi$ is
\begin{equation}\label{lim.final.qa.1}
 \lim_{w\to\infty}(\log w)\left(\frac{p-2}{n+2}- \phi(w)\right)=b\frac{2(p+n)}{(n+2)^2}=b\beta_\star.
\end{equation}
Thus we see that the behavior of the generalized Bayes shrinkage function in the cases
of known and unknown scale for the related classes of priors
are in very close correspondence. Additionally admissibility/inadmissibility in the known scale
case corresponds exactly with quasi-admissibility/inadmissibility in the unknown scale case.
We conjecture, for this class of priors in the unknown scale case, that
quasi-admissibility/inadmissibility in fact corresponds to admissibility/inadmissibility.

\end{remark}

 \section{Concluding remark}
\label{sec:cr}
We have studied quasi-admissible and quasi-inadmissible Stein-type shrinkage estimators
in the problem of estimating the mean vector of a $p$-variate Normal distribution when
the covariance matrix is an unknown multiple of the identity. We have established
sharp boundary of the form
\begin{equation}
 \phi_\star(w)=\frac{p-2}{n+2}-\frac{\beta_\star}{\log w}
\end{equation}
where $\beta_\star=2(p+2)/(n+2)^2$.
Roughly, estimators with shrinkage function $\phi(w)$ ultimately less than $\phi_\star(w)$
are quasi-inadmissible, while those which ultimately shrink more are quasi-admissible.
We have also found generalized prior distributions of the form $(1/\sigma^2)\times (1/\sigma^p)G(\|\theta\|/\sigma)$
for which he resulting generalized Bayes estimators are asymptotically of the form
\begin{equation*}
 \left\{1-\left(\frac{p-2}{n+2}-\frac{b\beta_\star}{\log W}\right)\frac{1}{W}\right\}X
\end{equation*}
for any $b>0$, thus establishing a boundary behavior for this class of priors
between quasi-admissibility and quasi-inadmissibility.
We conjecture, for this class of priors, that
quasi-admissibility/inadmissibility in fact corresponds to admissibility/inadmissibility.


\section{Introduction}
\label{sec:intro}

Modern machine learning applications often involve high-dimensional datasets with large sample sizes, on which simple algorithms such as variants of online gradient descent are competitive with more complicated batch algorithms.
In addition to often being more practical computationally as compared to batch methods, online methods apply to a wider range of scenarios such as online prediction against individual sequences.
In the online non-stochastic (adversarial) setting, the regret of a learning algorithm is a more natural quantity to analyze than the cumulative loss (for the latter an adversary may as well emit a data sequence of pure noise).
Although vanilla online gradient descent \citep{zinkevich2003online} obtains provably optimal regret in a minimax sense \citep{abernethy2008optimal}, if a method can simultaneously admit better guarantees against easy data sequences (such as low-dimensional data or i.i.d.~stochastic data) while maintaining the fallback guarantee of minimax optimality against the nastiest data sequences, such a method is even better. Such methods fall within the recently-sculpted space of ``learning faster from easy data'' \citep{grunwald2013learning}.

AdaGrad \citep{duchi2011adaptive} embodies an adaptive gradient family of algorithms, the two most prominent cases being full-matrix AdaGrad and diagonal AdaGrad. Whereas the former is computationally intractable in high-dimensions, the latter is tractable with only linear complexity in the dimension of the data. Both methods admit sequence-dependent regret bounds that can be much better than the regret bounds of online gradient descent. However, the computational levity of the diagonal version is matched with a price: if the components of the gradient are highly correlated, the diagonal version may fail to adapt well while the full-matrix version continues to be adaptive.

In addition to its theoretical strengths, even the diagonal version of AdaGrad has exhibited strong empirical performance on real-world problems \citep{duchi2011adaptive}. Moreover, diagonal AdaGrad also has been incorporated into deep learning algorithms that currently achieve state-of-the-art learning performance on a number of difficult image classification tasks \citep{dean2012large}.

Even though it is theoretically superior in terms of the regret, the full-matrix version has not seen similar large-scale applications due to its high per-round complexity of $O(n^3)$. This raises the question: is it possible to design a method lying between the full-matrix and diagonal versions of AdaGrad whose per-round complexity is $O(n)$? In this work, we answer this question in the affirmative, up to log factors, by presenting CompAdaGrad.
This method replaces full-matrix AdaGrad's Bregman divergence regularization with the sum of a \emph{compressed} Bregman divergence operating in a low-dimensional subspace and a diagonal Bregman divergence operating in the orthogonal complement of that subspace. The precise form can be seen in \eqref{eqn:psi-rp-oc}. This method admits a theoretical guarantee that appears to be between full-matrix AdaGrad and diagonal AdaGrad, as shown in \cref{sec:regret}.

It is natural to ask why one would compress at the regularization level rather than simply compressing the data itself via a random projection from the very beginning. There are several reasons why compression at the regularization level makes more sense.
\begin{enumerate}
\item An initial compression of the data can reduce computational complexity but may lose information irretrievably, whereas by shifting the compression into the regularization component of the objective some part of the otherwise lost complementary information can still be exploited.
\item In an adversarial setting, if the learning algorithm commits to a fixed low-dimensional subspace, even an oblivious adversary can ensure that all the interesting action occurs in the complement of this subspace.
\item In some applications, one actually needs a predictor in the original space. This might be for interpretability, for instance.
\end{enumerate}

In the next section, we review AdaGrad and derive CompAdaGrad. In \cref{sec:regret} we present a regret bound for CompAdaGrad. We show how to compute updates for certain composite regularizers in \cref{sec:computations}. This section also contains a result that may be of independent interest: \cref{thm:1-sparse} establishes $O(n)$ complexity for computing the $n$-dimensional Walsh-Hadamard Transform of a 1-sparse vector.
In \cref{sec:experiments}, we present experimental results on several datasets. Finally, we conclude the paper.

\section{AdaGrad: The Full, the Diagonal, and the Compressed}
\label{sec:adagrad}

\paragraph{Notation.}
Throughout this paper, we let $\mathcal{X} = \mathbb{R}^n$ and leave extensions to convex subsets for future work.
For a strongly convex, differentiable function $\psi$, let $B_\psi$ be the Bregman divergence induced from $\psi$, defined as
\begin{align*}
B_\psi(x, y) = \psi(x) - \psi(y) - \langle \nabla \psi(y), x - y \rangle .
\end{align*}

For a sequence of subdifferentiable convex functions $f_1, \ldots, f_T$, let $g_1, \ldots, g_T$ be a corresponding sequence of subgradients, so that for each $t \in [T]$ we have $g_t \in \partial f_t$.

\paragraph{Online Convex Optimization game.}

We consider the following online learning protocol parameterized by a convex regularization function $\varphi: \mathcal{X} \rightarrow \mathbb{R}_+$.

Let Nature be an oblivious adversary; that is, before the game begins Nature selects its sequence of functions possibly with knowledge of Learner's (potentially randomized) strategy. The game then proceeds over a sequence of rounds:

\begin{quote}
For round $t = 1, 2, ...$
\begin{enumerate}[label=(\arabic*)]
\item Learner makes a prediction $x_t$ in action space $\mathcal{X}$.
\item Nature reveals a convex loss function $f_t : \mathcal{X} \rightarrow \mathbb{R}_+$ .
\item Learner suffers composite loss $f_t(x_t) + \varphi(x_t)$.
\end{enumerate}
\end{quote}

The goal is to find a strategy that minimizes the regret, defined as follows.
Let $f_1, \ldots, f_T$ be a sequence of functions chosen by an oblivious adversary.
Then the regret on this sequence is
\begin{align} \label{eqn:regret}
\mathcal{R}(f_1, \ldots, f_T)
= \sum_{t=1}^T \bigl( f_t(x_t) + \varphi(x_t) \bigr)
    - \inf_{x^* \in \mathcal{X}} \sum_{t=1}^T \bigl( f_t(x^*) + \varphi(x^*) \bigr) .
\end{align}

AdaGrad with composite mirror descent embodies one family of learning strategies for obtaining low regret. Although previously AdaGrad also has been presented with regularized dual averaging, in this work we restrict to composite mirror descent for simplicity. AdaGrad for composite mirror descent is described by the updates
\begin{align} \label{eqn:update}
x_{t+1} = \argmin_{x \in \mathcal{X}}
              \left\{
                  \eta \langle g_t, x \rangle + \eta \varphi(x) + B_{\psi_t}(x, x_t)
               \right\} ,
\end{align}
for some constant learning rate $\eta > 0$ and some adaptive choice of convex function $\psi_t$.

We define $G_t := \sum_{s=1}^t g_s g_s^T$ and use the notation $\|x\|_A^2 = x^T A x$.
In its full-matrix incarnation, AdaGrad uses the choice \citep{duchi2011adaptive}
\begin{align*}
\psi_t(x) = \frac{1}{2} \| x \|_{G_t^{1/2}}^2 .
\end{align*}
As shown by \cite{duchi2011adaptive} and reproduced in \cref{sec:regret} for convenience, full-matrix AdaGrad admits a strongly adaptive regret bound. Unfortunately, the update \eqref{eqn:update} for full-matrix AdaGrad is not tractable for large $n$ because it involves a matrix square root and solving an $n$-dimensional linear system, each of which costs time $O(n^3)$. In response to this issue, there is a diagonal version of AdaGrad that admits updates in time $O(n)$. For a square matrix $A$, let $\mathrm{diag}(A)$ be the diagonal matrix satisfying $\mathrm{diag}(A)_{ii} = A_{ii}$ for all $i$ and $\mathrm{diag}(A)_{ij} = 0$ for all $(i,j)$ such that $i \neq j$. The diagonal version of AdaGrad uses the choice \citep{duchi2011adaptive}
\begin{align*}
\psi_t(x) = \frac{1}{2} \| x \|_{\mathrm{diag}(G_t)^{1/2}}^2 .
\end{align*}
This method also admits theoretical guarantees (cf.~\citep{duchi2011adaptive} or \cref{sec:regret}). However, by design the diagonal method ignores the correlations between the components of the gradients.

\paragraph{CompAdaGrad.}
In this work, we introduce \emph{CompAdaGrad}, a method which combines the full-matrix approach in a subspace plus the diagonal approach in the complementary subspace.

The idea of CompAdaGrad starts by restricting full-matrix AdaGrad's Bregman divergence to a low dimensional subspace by way of a mapping $\Pi: \mathbb{R}^n \rightarrow \mathbb{R}^k$ for some $k \leq n$ (and typically $k \ll n$). Since $\Pi G_t \Pi^T = \sum_{s=1}^T \Pi g_s (\Pi g_s)^T$, this leads to the modified Bregman divergence term
\begin{align} \label{eqn:radagrad}
\frac{1}{2} \| \Pi (x - x_t) \|_{(\Pi G_t \Pi^T)^{1/2}}^2
= \frac{1}{2} \| x - x_t \|_{\Pi^T (\Pi G_t \Pi)^{1/2} \Pi}^2 ,
\end{align}
which was also used by \cite{krummenacher2014radagrad}.\footnote{Those authors select $\Pi$ to be an SRHT, described below, and we do this as well.}
Our first remark is that if $k = n$ and $\Pi$ is in the orthogonal group, then it is easy to see that setting $\psi_t$ as in the RHS of \eqref{eqn:radagrad} recovers full-matrix AdaGrad.

A deficiency of \eqref{eqn:radagrad} is that it ignores all of the action in the orthogonal complement of the image of $\Pi$; however, this action can be addressed naturally by directly considering the action in this complementary subspace. To this end, define $P$ to be the orthogonal projector corresponding to $\Pi$, defined as $P := \Pi^T (\Pi \Pi^T)^{-1} \Pi$. The corresponding complementary orthogonal projector is then $P^\perp := I - P$. Incorporating complementary regularization into the Bregman divergence with a diagonal approximation (since the dimension of the complementary subspace $n - k$ is presumably high) yields
\begin{align} \label{eqn:radagrad-complementary}
\frac{1}{2} \| \Pi (x - x_t) \|_{(\Pi G_t \Pi^T)^{1/2}}^2
+ \frac{\tau}{2} \| P^\perp (x - x_t) \|_{\mathrm{diag}(G^\perp_t)^{1/2}}^2 ,
\end{align}
where $G^\perp_t := \sum_{s=1}^t (P^\perp g_s) (P^\perp g_s)^T = P^\perp G_T P^\perp$ and $\tau \geq 0$ is a parameter.

Note that when $k = n$ and $\Pi$ is in the orthogonal group, \eqref{eqn:radagrad-complementary} still recovers full-matrix AdaGrad (since the complementary subspace is empty). Additionally, when $k = 0$ and $\tau = 1$, \eqref{eqn:radagrad-complementary} recovers diagonal AdaGrad since $P^\perp = I$. Modulo a small modification to ensure strong convexity of $\psi_t$, we have just derived the method that we call CompAdaGrad, which is characterized by setting $\psi_t$ as
\begin{align} \label{eqn:psi-rp-oc}
\psi_t(x) = \frac{1}{2} \| x \|_{A^{(r)}_t + \tau A^{(c)}_t}^2
\end{align}
\begin{flalign*}
\text{for} &&
A^{(r)}_t := \Pi^T \left( \Pi G_t \Pi^T  + \delta_r I \right)^{1/2} \Pi &&
\text{and} &&
A^{(c)}_t :=
P^\perp
 \left( \mathrm{diag} \left( G^\perp_t \right)^{1/2} + \delta_c I \right)
P^\perp . &&
\end{flalign*}

\paragraph{SRHT specialization.}
Throughout the rest of this paper, we take $\Pi$ to be a Subsampled Randomized Hadamard Transform (SRHT), defined as $\Pi: = \sqrt{\frac{n}{k}} R H \Sigma$ for a row selector $R \in \mathbb{R}^{k \times n}$, an $n$-dimensional (orthogonal) Walsh-Hadamard matrix $H$, and a diagonal Rademacher matrix $\Sigma \in \mathbb{R}^{n \times n}$ (whose diagonal entries are drawn i.i.d.~ as $+1$ and $-1$ with equal probability). The row selector $R$ is induced from the distribution of all cardinality-$k$ subsets of $n$ indices.
The idea behind \eqref{eqn:psi-rp-oc} is that if there is interesting action in a low-dimensional subspace, the SRHT (i.e.~the $A^{(r)}$-part of the regularization) potentially can capture this action \citep{tropp2011improved}.
However, as the SRHT can only capture action in a low-dimensional subspace, it is critical to also include \emph{complementary} regularization by way of the $A^{(c)}$ part of the regularization. This is especially true in the adversarial regime, where an adversary can pick up on Learner's fixed SRHT and make all of the interesting action happen in the complement. In this sense, CompAdaGrad can have the benefit of full-matrix AdaGrad under an oblivious, easy adversary who restricts most of the action in a low-dimensional subspace, while also maintaining guarantees against a harder adversary who uses their knowledge of $\mathrm{Im}(P)$.

The choices of $k$ and $\tau$ offer useful degrees of freedom for massive datasets. As $k$ increases toward $n$, the method more closely resembles full-matrix AdaGrad, and the regret bound generally becomes stronger as $k$ increases, while the computational complexity for various methods also becomes larger with larger $k$. The parameter $\tau$ on the other hand allows one to modulate how much emphasis to place on the complementary subspace, with this action increasingly ignored as $\tau$ decreases to 0. The next two sections help provide an understanding of the trade-off between regret and computation respectively.

\section{Regret Bounds}
\label{sec:regret}

Diagonal AdaGrad inherently is unable to adapt to correlations in a data sequence, as can be seen by Duchi et al.'s regret bound for this method (cf.~Theorem 5 of \cite{duchi2011adaptive}):
\begin{theorem} \label{thm:diagonal-regret}
\begin{align} \label{eqn:diagonal-regret}
\mathcal{R}(f_1, \ldots, f_T)
\leq \frac{1}{2 \eta} \max_{t \in [T]} \|x^* - x_t\|_\infty^2
                     \sum_{j=1}^n \left( \sum_{t=1}^T g_{t,j}^2 \right)^{1/2}
                     + \eta \sum_{j=1}^n \left( \sum_{t=1}^T g_{t,j}^2 \right)^{1/2}
\end{align}
\end{theorem}
Although diagonal AdaGrad can perform well on high-dimensional sparse data, the above summations over the $n$ dimensions are symptomatic of its inability to adapt to highly correlated dimensions.

In contrast, full-matrix AdaGrad admits the regret bound (cf.~Theorem 7 of \cite{duchi2011adaptive})
\begin{theorem} \label{thm:full-regret}
\begin{align} \label{eqn:full-regret}
\mathcal{R}(f_1, \ldots, f_T)
\leq \frac{\delta}{\eta} \|x^*\|_2^2
                     + \frac{1}{2 \eta} \max_{t \in [T]} \|x^* - x_t\|_2^2 \tr(G_T^{1/2})
                     + \eta \tr(G_T^{1/2}) .
\end{align}
\end{theorem}
The above result, depending primarily on the trace of the covariance of the gradients, can exploit high correlations among the dimensions of points in a data sequence. We will show that, in certain situations of interest, CompAdaGrad can come close to the above regret guarantee at a small fraction of the computational complexity.

By design, it is quite straightforward to work out a regret bound for CompAdaGrad by leveraging the existing analysis of both full-matrix and diagonal AdaGrad from \cite{duchi2011adaptive}. Although in the context of concentration inequalities it is important to use the SRHT scaled as $\Pi = \sqrt{\frac{n}{k}} R H \Sigma$, in our analysis we instead analyze a variant of CompAdaGrad that is defined by replacing $\Pi$ with the unscaled $\tilde{\Pi} := R H \Sigma$ (which notably satisfies $\tilde{\Pi} \tilde{\Pi}^T = I$). The following regret bound is for this variant of CompAdaGrad. Also, we set $\delta_r = \delta_c = \delta$ because it simplifies the presentation.
\begin{theorem} \label{thm:compressed-regret}
CompAdaGrad with learning rate $\eta > 0$ and $\delta > 0$ satisfies
\begin{align} \label{eqn:compressed-regret}
\mathcal{R}(f_1, \ldots, f_T)
\leq
&\frac{\delta}{2 \eta} \|x^* - x_1\|_2^2 \nonumber \\
&\frac{1}{2 \eta}
        \left(
            \max_{t \in [T]} \| P (x^* - x_t) \|_2^2
            \tr \left( (\tilde{\Pi} G_T \tilde{\Pi}^T)^{1/2} \right)
            + \max_{t \in [T]} \| P^\perp (x^* - x_t) \|_\infty^2 \| Z^\perp \|_{2,1}
        \right) \nonumber \\
       &+ \eta \left(
               \tr \left( (\tilde{\Pi} G_T \tilde{\Pi}^T)^{1/2} \right)
               + \| Z^\perp \|_{2,1}
           \right) ,
\end{align}
where $Z^\perp = (P^\perp g_1 \, \ldots \, P^\perp g_T)^T$ and $A \mapsto \| A \|_{2,1}$ is the sum of the $\ell_2$ norms of the columns of $A$.
\end{theorem}
The proof can be found in \cref*{appendix:regret}.

Note that in the case where the data occupies a low-dimensional subspace and the SRHT preserves the action within this subspace, the above bound is similar to the regret bound for full-matrix AdaGrad \eqref{eqn:full-regret}. However, even when the data sequence is not so easy, \cref{thm:compressed-regret} still offers a fall-back guarantee based on the action in the orthogonal complement. In the event that the dimensions of the gradients are uncorrelated and the gradients do not occupy a low-dimensional subspace, the guarantee for diagonal AdaGrad provided by \cref{thm:diagonal-regret} could be better than the guarantee for CompAdaGrad provided by \cref{thm:compressed-regret}.

\section{Computations}
\label{sec:computations}

In this section, we show two important composite regularizers for which it is possible to compute the updates steps for CompAdaGrad: the squared $\ell_2$ regularizer and the $\ell_1$ regularizer.
Some of the results below rely upon the following conjecture which we hope to affirm in the long version.
\begin{conjecture} \label{conjecture:k-sparse}
The $n$-dimensional Walsh-Hadamard Transform of an $r$-sparse vector can be computed in time $O(n \log r)$.
\end{conjecture}
For any results that rely on the conjecture for some variable (e.g.~$r$) indicating the sparsity level, we present the results instead with a \emph{primed} version of the variable (e.g.~$r'$) with the understanding that $r' = r$ if the conjecture is true and $r' = n$ otherwise.

Regardless of the veracity of \cref{conjecture:k-sparse}, we do however prove the following weaker result.
\begin{theorem} \label{thm:1-sparse}
The $n$-dimensional Walsh-Hadamard Transform of a $1$-sparse vector can be computed in time $O(n)$.
\end{theorem}
The proof is constructive, and the algorithm and its analysis can be found in \cref*{appendix:1-sparse-proof}.

\subsection{CompAdaGrad with the squared $\ell_2$ composite regularizer}

An update for squared $\ell_2$ composite mirror descent with the compressed Bregman divergence can be written in the form
\begin{align*}
\min_x \eta \langle g_t, x \rangle
            + \frac{1}{2} \left\| \Pi (x- x_t) \right\|_{K_t}^2
            + \frac{\tau}{2} \left\| P^\perp (x - x_t) \right\|_{D_t}^2
            + \frac{\eta \lambda}{2} \|x\|_2^2 ,
\end{align*}
where $K_t = (\Pi G_t \Pi^T + \delta_r I)^{1/2}$ and $D_t = \mathrm{diag}(G^\perp_t)^{1/2} + \delta_c I$.

To simplify, we define $g := \eta g_t$, $K := K_t$, and $D = \tau D_t$, and we replace $\eta \lambda$; the above is then
\begin{align} \label{eqn:squared-ell2-simplified}
\min_x \langle g, x \rangle
            + \frac{1}{2} \left\| \Pi (x- x_t) \right\|_K^2
            + \frac{1}{2} \left\| P^\perp (x - x_t) \right\|_D^2
            + \frac{\lambda}{2} \|x\|_2^2 .
\end{align}
Since $P$ is the orthogonal projector corresponding to the low-dimensional mapping $\Pi$, the above can be written equivalently as
\begin{align*}
\min_x \left\{
  \textcolor{dark-blue}{
    \langle P g, x \rangle
    + \frac{1}{2} \left\| P (x- x_t) \right\|_{\Pi^T K \Pi}^2
    + \frac{\lambda}{2} \|P x\|_2^2
  }
  + \textcolor{dark-red}{
    \langle P^\perp g, x \rangle
    + \frac{1}{2} \left\| P^\perp (x - x_t) \right\|_D^2
    + \frac{\lambda}{2} \|P^\perp x\|_2^2
  }
\right\} .
\end{align*}

From the above final rewriting, it is clear that the optimization decouples over the two subspaces, with the minimizer of the above problem $x^*$ being equal to the sum of
\begin{align} \label{eqn:problem-parallel}
&&& x^*_\parallel :=
  \argmin_{x \in \mathrm{Im}(P)} \left\{
    \textcolor{dark-blue}{
       \langle P g, x \rangle
      + \frac{1}{2} \left\| x - P x_t \right\|_{\Pi^T K \Pi}^2
      + \frac{\lambda}{2} \|x\|_2^2
    }
\right\} && \\
\text{and} &&& x^*_\perp :=
\argmin_{x \in \mathrm{Im}(P^\perp)} \left\{
  \textcolor{dark-red}{
    \langle P^\perp g, x \rangle
    + \frac{1}{2} \left\| x - P^\perp x_t \right\|_D^2
    + \frac{\lambda}{2} \|x\|_2^2
   }
    \right\} . && \label{eqn:problem-perp}
\end{align}

As shown in \cref*{appendix:computational}, the solution to \eqref{eqn:problem-parallel} is
\begin{align}
x^*_\parallel = (\Pi^T K \Pi + \lambda I)^{-1} (\Pi^T K \Pi x_t - P g) . \label{eqn:unconstrained-constrained}
\end{align}
This can be computed efficiently as follows; the proof can be found in \cref*{appendix:computational}.
\begin{lemma} \label{lemma:computation-squared-ell2-parallel}
$x^*_\parallel $ can be computed as
\begin{align} \label{eqn:computation-squared-ell2-parallel}
x^*_\parallel
&= \Pi^T \left( \frac{n}{k} K + \lambda I \right)^{-1} \left( K \Pi x_t - \frac{k}{n} \Pi g \right)
\end{align}
in time $O(n \log k' + k^3)$.
\end{lemma}

To solve \eqref{eqn:problem-perp}, we take the dual. First, observe that \eqref{eqn:problem-perp} can be rewritten as
\begin{align*}
\begin{aligned}
&\min_x
& & \langle P^\perp g, x \rangle
       + \frac{1}{2} \left\| x - P^\perp x_t \right\|_D^2
       + \frac{\lambda}{2} \|x\|_2^2 \\
&\,\,\text{s.t.}
& & \Pi x = \mathbf{0} .
\end{aligned}
\end{align*}

The Lagrangian is
\begin{align*}
\mathcal{L}(x, \nu) =
  \langle P^\perp g, x \rangle
  + \frac{1}{2} \left\| x - P^\perp x_t \right\|_D^2
  + \frac{\lambda}{2} \|x\|_2^2
  + \langle \nu, \Pi x \rangle .
\end{align*}

\begin{lemma} \label{lemma:computation-squared-ell2-perp}
Problem \eqref{eqn:problem-perp} can be solved as
\begin{align} \label{eqn:x-in-terms-of-nu}
x^*_\perp = (D + \lambda I)^{-1} (-P^\perp g + D P^\perp x_t - \Pi^T \nu^*)
\end{align}
with
\begin{align} \label{eqn:computation-for-nu}
\nu^* =
\left( \Pi B \Pi^T \right)^{-1}
\Pi B (-P^\perp g + D P^\perp x_t - D B y - \lambda B y + y) ,
\end{align}
for $B = (D + \lambda I)^{-1}$ and $y = -P^\perp g + D P^\perp x_t$.
These computations can be completed in time $O(n k \log k)$.
\end{lemma}
The proof can be found in \cref*{appendix:computational}.

\subsection{CompAdaGrad with the $\ell_1$ composite regularizer}

We first sketch the LARS algorithm \citep{efron2004least}. We then show how the computations can be done efficiently for CompAdaGrad.

\paragraph{Re-expression as a LASSO problem.}

Using the same simplifications as in \eqref{eqn:squared-ell2-simplified}, the update step is
\begin{align*}
\argmin_x \langle g, x \rangle
            + \frac{1}{2} \left\| x- x_t \right\|_A^2
            + \lambda \|x\|_1 .
\end{align*}
for the symmetric matrix $A = \Pi^T K \Pi + P^\perp D P^\perp$. The above can be re-expressed as
\begin{align} \label{eqn:l1-update}
  x_{t+1} = \argmin_x \left\{
                                      \langle u, x \rangle
                                      + \frac{1}{2} \langle x, A x \rangle
                                      + \lambda \|x\|_1
                                  \right\} ,
\end{align}
for $u = g - A x_t$. Note that $A$ is positive definite whenever $\delta_r, \delta_c > 0$.

The problem \eqref{eqn:l1-update} can be expressed as a LASSO problem, as the optimal objective value is equal to
\begin{align*}
\min_x \frac{1}{2} \|x\|_A^2 - \langle -A^{-1} u, x \rangle_A + \frac{1}{2} \left\| -A^{-1} u \right\|_A^2 + \lambda \|x\|_1
&\equiv
\min_x \frac{1}{2} \left\| x - (-A^{-1} u) \right\|_A^2 + \lambda \|x\|_1 \\
&\equiv
\min_x \frac{1}{2} \left\| A^{1/2} x - (-A^{-1/2} u) \right\|^2 + \lambda \|x\|_1 .
\end{align*}

\paragraph{Efficient computations}

The standard LARS algorithm takes as input a matrix of covariates $X$ and targets $y$. In our case, $X = A^{1/2}$ is too expensive to compute since $A \in \mathbb{R}^{n \times n}$, while $y = A^{-1/2} u$ is similarly too expensive to compute. Therefore, we cannot just run the usual LARS algorithm but instead have to find efficient ways to perform some of the algorithm's internal computations.

The two computations for LARS that need to be handled differently are the computation of
(a) the correlation of the covariates with the current residual;
and
(b) the entries of the Gram matrix for doing Cholesky Insert operations.

For both of the above, a certain $k$-by-$k$ matrix $Q$ arises in the computations which can be pre-computed and re-used for the entire run of LARS. $Q$ is defined as
$\frac{n}{k} K + R H \Sigma D \Sigma H R^T$, or equivalently as $\frac{n}{k} K + R H D H R^T$. The second term can computed in time $O(n k \log k)$, in parallel over columns, precisely like $\Pi B \Pi^T$ from \eqref{eqn:computation-for-nu}, as explained in the proof of \cref{lemma:computation-squared-ell2-perp}.

Given a current hypothesis $\beta$, the correlation of the covariates with the residual can be computed as $X^T(y - X \beta) = A^{1/2} A^{-1/2} u - A^{1/2} A^{1/2} \beta = u - A \beta$. Thus, a matrix vector multiplication is the main work. This operation can be completed in time $O(n \log k' + k^2)$, as shown in \cref*{appendix:matrix-vector}.

When performing a Cholesky Insert (increasing the active set $\Lambda$ by one), we need to compute part of a column of the Gram matrix. If the new dimension is $i$, we need to compute $G_{ii}$ and $G_{ij}$ for each $j \in \Lambda$. Since the Gram matrix $G = A$, we just need to compute $e_i^T A e_j$ for each $j \in \Lambda$.

To compute $G_{ii}$, observe that
\begin{align} \label{eqn:Gram-ii}
G_{ii} = A_{ii} = d_i \langle R H e_i, Q R H e_i \rangle - 2 d_i \frac{k}{n} ,
\end{align}
To compute \eqref{eqn:Gram-ii}, we first compute $R H e_i$ in $O(k \log n)$ and store this for use in future rounds. The additional cost for computing $Q R H e_i$ is then $O(k^2)$, and we re-use this result in \eqref{eqn:Gram-column}. So the cost is $O(k \log n + k^2)$.
To compute $G_{i,\Lambda}$, we compute for each $j \in \Lambda$
\begin{align} \label{eqn:Gram-column}
G_{ij} = \sigma_i \sigma_j \left( \langle R H e_j, Q R H e_i \rangle - (d_i + d_j) \langle R H e_j, R H e_i \rangle \right) .
\end{align}
Computing \eqref{eqn:Gram-column} can be done in $O(|\Lambda| k)$ since we have already stored $R H e_j$ for all $j \in \Lambda$ and already computed $Q R H e_i$.

\section{Experiments}
\label{sec:experiments}

For all the experiments we used the squared $\ell_2$ regularizer. For CompAdaGrad, we always set $\delta_r$ and $\delta_c$ to be equal and hereafter refer to them both as $\delta$.
When reporting test error for MNIST and Reuters RCV1, we select the values of $\delta$, $\eta$, $\lambda$, and $\tau$ that attained the lowest online zero-one loss.

\paragraph{MNIST 4 vs 9.}
This dataset was constructed by selecting 400 random prototypes (200 from each class) from the training set and using as features a Gaussian kernel computation of each data point with each of the 400 prototypes. We train on the entire training set for MNIST 4 and 9 and test on the test set for MNIST 4 and 9.
The results in \cref{fig:mnist} indicate that CompAdaGrad well-outperforms the diagonal method, and as the SRHT dimension increases from 25 to 256, the gains over the diagonal method are all the larger. It is worth mentioning that the CompAdaGrad improves over diagonal AdaGrad even for $k = 25$. We suspect that that CompAdaGrad surpasses diagonal AdaGrad on this task because the features are dense and highly correlated.

\paragraph{Reuters RCV1.}
We used Reuters RCV1-v2 \citep{lewis2004rcv1}. After stopping, stemming, and retaining only those unigrams and bigrams that occurred at least twice in the corpus, the data consisted of 800,000 documents in 1,889,478 dimensions. We computed gradients in mini-batches of 160 points and performed a hypothesis update from each resulting averaged gradient. We generated 4 random permutations of the data. For each permutation, we train on the first 75\% and test on the last 25\%. Each of the 4 runs is an experiment.
We report the average test error over the 4 experiments. This test error need not correspond to a single $(\eta, \lambda,\delta,(\tau))$ configuration.
As shown in \cref{table:reuters}, CompAdaGrad with $k=512$ does not improve upon the diagonal version. We suspect the lack of improvement is because the data is very sparse and does not admit a low-dimensional linear subspace capturing a large amount of the action in the data. Therefore, it does not help much to pick up on correlations.

\paragraph{Foxes and Wolves.}
We grabbed a visual-bag-of-words representation for 829 images of kit foxes (\emph{Vulpes macrotis}) and 1,156 images of red wolves (\emph{Canis rufus}) from ImageNet \citep{russakovsky2015imagenet} (originally from Flickr). The original features are 1000 visual words which represent 1000 clusters of a random subset of 10 million SIFT features.
We computed visual word bigrams by selecting\footnote{We selected 0.002 based on the quantiles of the distance distribution of all pairs of visual words in images.} a radius of 0.002 for visual words to be considered collocated (based on $(x,y)$ coordinates for the visual words); this led to 238,822 bigrams which occurred at least once in the 1,985 images. The final dataset for learning consisted of the unigram counts and bigram counts, constituting 239,822 features.

We trained on the first 75\% of a random permutation of the data and tested on the remainder.
Because the number of instances in this dataset is small relative to the dimension, the online zero-one loss is not a sensible rule for parameter tuning as a relatively large number hypotheses contributing to the online loss are from early rounds where not much learning has taken place. We therefore tuned the parameters based on the online training zero-one loss of the final hypothesis of the online learning algorithm. The results presented in \cref{table:foxes-wolves} are inconclusive. Both diagonal AdaGrad and CompAdaGrad with $k=256$ obtain the same zero-one test risk. This either could be due to the hardness of the task, the small sample size, or the fact that the features are still quite sparse due to the visual bag-of-words bigram representation.

\begin{figure}
\begin{floatrow}
\capbtabbox{
\begin{tabular}{c c}
\toprule
Method                   & Zero-one test risk \\
\midrule
Diagonal            &  0.056084 \\
CompAdaGrad-512  & 0.056155 \\
\bottomrule
\end{tabular}
}{
  \caption{\label{table:reuters} Reuters RCV1 test error.}
}
\capbtabbox{
\begin{tabular}{c c}
\toprule
Method      & Zero-one test risk \\
\midrule
Diagonal   &  0.2696 \\
CompAdaGrad-256 & 0.2696 \\
\bottomrule
\end{tabular}
}{
  \caption{\label{table:foxes-wolves}  Foxes and Wolves test error}
}
\end{floatrow}
\end{figure}

\begin{figure}
\begin{floatrow}
\ffigbox{
\includegraphics[width=75mm]{mnist_k_vs_test_error.pdf}
}{
  \caption{\label{fig:mnist} MNIST 4 vs 9 test error.}
}
\ffigbox{
\includegraphics[width=40mm]{kit_fox.jpg}
\includegraphics[width=40mm]{red_wolf.jpg}
}{
  \caption{A kit fox and a red wolf}
}
\end{floatrow}
\end{figure}

\section{Conclusion}
\label{sec:discussion}

We have introduced a new adaptive gradient method, CompAdaGrad, which can obtain regret bounds competitive with full-matrix AdaGrad under easier adversaries that put most of the interesting action in a low-dimensional subspace, while maintaining a fallback regret bound in the case that the adversary is not easy. CompAdaGrad also admits efficient updates for certain choices of the composite regularizer. In preliminary experiments, we demonstrate that CompAdaGrad can sometimes obtain lower risk as compared to diagonal AdaGrad, but this outcome appears to be tied to the density and more importantly the correlation of the features. In the future, we intend to apply CompAdaGrad to massive datasets with dense, highly correlated features. Along these lines, one promising application for future work is to incorporate CompAdaGrad into the training of deep learning methods.

\subsubsection*{Acknowledgments}

This work is supported in part by the Australian Research Council Discovery Project DP0987773 and NICTA, which is funded by the Australian Government through the Department of Communications and the Australian Research Council through the ICT Centre of Excellence Program.

\newpage

\renewcommand{\refname}{{\normalsize References}}