Randomized Structural Sparsity based Support Identification with Applications to Locating Activated or Discriminative Brain Areas: A Multi-center Reproducibility Study
Yilun~Wang
%~\IEEEmembership{Member
~IEEE
....lastname \thanks{...
In this paper, we focus on how to locate the relevant or discriminative brain regions related with external stimulus or certain mental decease, which is also called support identification, based on the neuroimaging data. The main difficulty lies in the extremely high dimensional voxel space and relatively few training samples, easily resulting in an unstable brain region discovery (or called feature selection in context of pattern recognition). When the training samples are from different centers and have between-center variations, it will be even harder to obtain a reliable and consistent result. %for many traditional multivarate feature selection methods. %In order to deal with these difficulties, many efforts have been performed and we will mainly focus on designing stable feature selection algorithms. Corresponding, we revisit our recently proposed algorithm based on stability selection and % % %stability selection for discriminative features have recently received great attentions, due to its %finite sample control for certain error rates of false discoveries and %transparent principle for choosing a proper amount of regularization. %However, it fails to consider some prior or estimated structural information of these features. Correspondingly, we propose a %algorithm which is able to incorporate the prior % structural sparsity. % is able to achieve stable and reliable support identification. %While its advantages over the common stability selection is the better control of false negatives, yet keeping the control of false positives. %Compared with the common structural sparsity models %, the prior or estimated structural information is only for a reference and not necessarily very accurate, which is important in practice. It is applied to the multi-center MRI data analysis for the first time. A consistent and stable result is achieved across different centers despite the between-center data variation while many other state-of-the-art methods such as two sample t-test fail. Moreover, we have empirically showed that the performance of this algorithm is robust and insensitive to several of its key parameters. In addition, the support identification results on both functional MRI and structural MRI are interpretable and can be the potential biomarkers. %In this paper, we focus on the functional magnetic resonance imaging data of the brain, and we present a specific implementation based on our proposed %I. % %voxels from functional magnetic resonance imaging data have been gaining great attention recently, because the widely used univariate methods such as the %two-sample t-test fail to find multi-voxel patterns. %, and remain an unsolved problem, because more attention is paid to the prediction accuracy instead. %In this study, we focus on sparse pattern recognition methods %based identification of the involved brain regions in %for very high dimensional functional %MRI (fMRI) data associated with either a cognitive task or disease diagnosis. %However, the inferred discriminative voxels of most existing methods %are instantiated and optimized for the prediction accuracy, but their inferred discriminative voxels %either have significant false negatives or false positives, %However, most of them fail to give accurate and stable results in the cases of very high dimensional voxel space. In order to partially overcome this problem, we proposed a variant of stability selection based on % %``constrained block subsampling" to %and multi-stage refineness aiming to %it is able to reduce the falsely selected voxels, aiming to % reliably and stably recover the ground truth of the discriminative voxels efficiently. %Numerical experiments %We % demonstrate that our method can achieve better accuracy and reproducibility than many alternatives based on the multi-center autism brain imaging data across different modals. %Specifically, compared with the original stability selection, this variant
bare_jrnl.tex

%% bare_jrnl.tex
%% V1.4
%% 2012/12/27
%% by Michael Shell
%% see http://www.michaelshell.org/
%% for current contact information.
%%
%% This is a skeleton file demonstrating the use of IEEEtran.cls
%% (requires IEEEtran.cls version 1.8 or later) with an IEEE journal paper.
%%
%% Support sites:
%% http://www.michaelshell.org/tex/ieeetran/
%% http://www.ctan.org/tex-archive/macros/latex/contrib/IEEEtran/
%% and
%% http://www.ieee.org/



% *** Authors should verify (and, if needed, correct) their LaTeX system  ***
% *** with the testflow diagnostic prior to trusting their LaTeX platform ***
% *** with production work. IEEE's font choices can trigger bugs that do  ***
% *** not appear when using other class files.                            ***
% The testflow support page is at:
% http://www.michaelshell.org/tex/testflow/


%%*************************************************************************
%% Legal Notice:
%% This code is offered as-is without any warranty either expressed or
%% implied; without even the implied warranty of MERCHANTABILITY or
%% FITNESS FOR A PARTICULAR PURPOSE!
%% User assumes all risk.
%% In no event shall IEEE or any contributor to this code be liable for
%% any damages or losses, including, but not limited to, incidental,
%% consequential, or any other damages, resulting from the use or misuse
%% of any information contained here.
%%
%% All comments are the opinions of their respective authors and are not
%% necessarily endorsed by the IEEE.
%%
%% This work is distributed under the LaTeX Project Public License (LPPL)
%% ( http://www.latex-project.org/ ) version 1.3, and may be freely used,
%% distributed and modified. A copy of the LPPL, version 1.3, is included
%% in the base LaTeX documentation of all distributions of LaTeX released
%% 2003/12/01 or later.
%% Retain all contribution notices and credits.
%% ** Modified files should be clearly indicated as such, including  **
%% ** renaming them and changing author support contact information. **
%%
%% File list of work: IEEEtran.cls, IEEEtran_HOWTO.pdf, bare_adv.tex,
%%                    bare_conf.tex, bare_jrnl.tex, bare_jrnl_compsoc.tex,
%%                    bare_jrnl_transmag.tex
%%*************************************************************************

% Note that the a4paper option is mainly intended so that authors in
% countries using A4 can easily print to A4 and see how their papers will
% look in print - the typesetting of the document will not typically be
% affected with changes in paper size (but the bottom and side margins will).
% Use the testflow package mentioned above to verify correct handling of
% both paper sizes by the user's LaTeX system.
%
% Also note that the "draftcls" or "draftclsnofoot", not "draft", option
% should be used if it is desired that the figures are to be displayed in
% draft mode.
%

\documentclass[11pt,draftcls,onecolumn] {IEEEtran}
%\documentclass[journal]{IEEEtran}
%
% If IEEEtran.cls has not been installed into the LaTeX system files,
% manually specify the path to it like:
% \documentclass[journal]{../sty/IEEEtran}





% Some very useful LaTeX packages include:
% (uncomment the ones you want to load)


% *** MISC UTILITY PACKAGES ***
%
%\usepackage{ifpdf}
% Heiko Oberdiek's ifpdf.sty is very useful if you need conditional
% compilation based on whether the output is pdf or dvi.
% usage:
% \ifpdf
%   % pdf code
% \else
%   % dvi code
% \fi
% The latest version of ifpdf.sty can be obtained from:
% http://www.ctan.org/tex-archive/macros/latex/contrib/oberdiek/
% Also, note that IEEEtran.cls V1.7 and later provides a builtin
% \ifCLASSINFOpdf conditional that works the same way.
% When switching from latex to pdflatex and vice-versa, the compiler may
% have to be run twice to clear warning/error messages.






% *** CITATION PACKAGES ***
%
%\usepackage{cite}
% cite.sty was written by Donald Arseneau
% V1.6 and later of IEEEtran pre-defines the format of the cite.sty package
% \cite{} output to follow that of IEEE. Loading the cite package will
% result in citation numbers being automatically sorted and properly
% "compressed/ranged". e.g., [1], [9], [2], [7], [5], [6] without using
% cite.sty will become [1], [2], [5]--[7], [9] using cite.sty. cite.sty's
% \cite will automatically add leading space, if needed. Use cite.sty's
% noadjust option (cite.sty V3.8 and later) if you want to turn this off
% such as if a citation ever needs to be enclosed in parenthesis.
% cite.sty is already installed on most LaTeX systems. Be sure and use
% version 4.0 (2003-05-27) and later if using hyperref.sty. cite.sty does
% not currently provide for hyperlinked citations.
% The latest version can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/cite/
% The documentation is contained in the cite.sty file itself.






% *** GRAPHICS RELATED PACKAGES ***
%
\ifCLASSINFOpdf
  % \usepackage[pdftex]{graphicx}
  % declare the path(s) where your graphic files are
  % \graphicspath{{../pdf/}{../jpeg/}}
  % and their extensions so you won't have to specify these with
  % every instance of \includegraphics
  % \DeclareGraphicsExtensions{.pdf,.jpeg,.png}
\else
  % or other class option (dvipsone, dvipdf, if not using dvips). graphicx
  % will default to the driver specified in the system graphics.cfg if no
  % driver is specified.
  % \usepackage[dvips]{graphicx}
  % declare the path(s) where your graphic files are
  % \graphicspath{{../eps/}}
  % and their extensions so you won't have to specify these with
  % every instance of \includegraphics
  % \DeclareGraphicsExtensions{.eps}
\fi
% graphicx was written by David Carlisle and Sebastian Rahtz. It is
% required if you want graphics, photos, etc. graphicx.sty is already
% installed on most LaTeX systems. The latest version and documentation
% can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/required/graphics/
% Another good source of documentation is "Using Imported Graphics in
% LaTeX2e" by Keith Reckdahl which can be found at:
% http://www.ctan.org/tex-archive/info/epslatex/
%
% latex, and pdflatex in dvi mode, support graphics in encapsulated
% postscript (.eps) format. pdflatex in pdf mode supports graphics
% in .pdf, .jpeg, .png and .mps (metapost) formats. Users should ensure
% that all non-photo figures use a vector format (.eps, .pdf, .mps) and
% not a bitmapped formats (.jpeg, .png). IEEE frowns on bitmapped formats
% which can result in "jaggedy"/blurry rendering of lines and letters as
% well as large increases in file sizes.
%
% You can find documentation about the pdfTeX application at:
% http://www.tug.org/applications/pdftex





% *** MATH PACKAGES ***
%
%\usepackage[cmex10]{amsmath}
% A popular package from the American Mathematical Society that provides
% many useful and powerful commands for dealing with mathematics. If using
% it, be sure to load this package with the cmex10 option to ensure that
% only type 1 fonts will utilized at all point sizes. Without this option,
% it is possible that some math symbols, particularly those within
% footnotes, will be rendered in bitmap form which will result in a
% document that can not be IEEE Xplore compliant!
%
% Also, note that the amsmath package sets \interdisplaylinepenalty to 10000
% thus preventing page breaks from occurring within multiline equations. Use:
%\interdisplaylinepenalty=2500
% after loading amsmath to restore such page breaks as IEEEtran.cls normally
% does. amsmath.sty is already installed on most LaTeX systems. The latest
% version and documentation can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/required/amslatex/math/





% *** SPECIALIZED LIST PACKAGES ***
%
%\usepackage{algorithmic}
% algorithmic.sty was written by Peter Williams and Rogerio Brito.
% This package provides an algorithmic environment fo describing algorithms.
% You can use the algorithmic environment in-text or within a figure
% environment to provide for a floating algorithm. Do NOT use the algorithm
% floating environment provided by algorithm.sty (by the same authors) or
% algorithm2e.sty (by Christophe Fiorio) as IEEE does not use dedicated
% algorithm float types and packages that provide these will not provide
% correct IEEE style captions. The latest version and documentation of
% algorithmic.sty can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/algorithms/
% There is also a support site at:
% http://algorithms.berlios.de/index.html
% Also of interest may be the (relatively newer and more customizable)
% algorithmicx.sty package by Szasz Janos:
% http://www.ctan.org/tex-archive/macros/latex/contrib/algorithmicx/




% *** ALIGNMENT PACKAGES ***
%
%\usepackage{array}
% Frank Mittelbach's and David Carlisle's array.sty patches and improves
% the standard LaTeX2e array and tabular environments to provide better
% appearance and additional user controls. As the default LaTeX2e table
% generation code is lacking to the point of almost being broken with
% respect to the quality of the end results, all users are strongly
% advised to use an enhanced (at the very least that provided by array.sty)
% set of table tools. array.sty is already installed on most systems. The
% latest version and documentation can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/required/tools/


% IEEEtran contains the IEEEeqnarray family of commands that can be used to
% generate multiline equations as well as matrices, tables, etc., of high
% quality.




% *** SUBFIGURE PACKAGES ***
%\ifCLASSOPTIONcompsoc
%  \usepackage[caption=false,font=normalsize,labelfont=sf,textfont=sf]{subfig}
%\else
%  \usepackage[caption=false,font=footnotesize]{subfig}
%\fi
% subfig.sty, written by Steven Douglas Cochran, is the modern replacement
% for subfigure.sty, the latter of which is no longer maintained and is
% incompatible with some LaTeX packages including fixltx2e. However,
% subfig.sty requires and automatically loads Axel Sommerfeldt's caption.sty
% which will override IEEEtran.cls' handling of captions and this will result
% in non-IEEE style figure/table captions. To prevent this problem, be sure
% and invoke subfig.sty's "caption=false" package option (available since
% subfig.sty version 1.3, 2005/06/28) as this is will preserve IEEEtran.cls
% handling of captions.
% Note that the Computer Society format requires a larger sans serif font
% than the serif footnote size font used in traditional IEEE formatting
% and thus the need to invoke different subfig.sty package options depending
% on whether compsoc mode has been enabled.
%
% The latest version and documentation of subfig.sty can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/subfig/




% *** FLOAT PACKAGES ***
%
%\usepackage{fixltx2e}
% fixltx2e, the successor to the earlier fix2col.sty, was written by
% Frank Mittelbach and David Carlisle. This package corrects a few problems
% in the LaTeX2e kernel, the most notable of which is that in current
% LaTeX2e releases, the ordering of single and double column floats is not
% guaranteed to be preserved. Thus, an unpatched LaTeX2e can allow a
% single column figure to be placed prior to an earlier double column
% figure. The latest version and documentation can be found at:
% http://www.ctan.org/tex-archive/macros/latex/base/


%\usepackage{stfloats}
% stfloats.sty was written by Sigitas Tolusis. This package gives LaTeX2e
% the ability to do double column floats at the bottom of the page as well
% as the top. (e.g., "\begin{figure*}[!b]" is not normally possible in
% LaTeX2e). It also provides a command:
%\fnbelowfloat
% to enable the placement of footnotes below bottom floats (the standard
% LaTeX2e kernel puts them above bottom floats). This is an invasive package
% which rewrites many portions of the LaTeX2e float routines. It may not work
% with other packages that modify the LaTeX2e float routines. The latest
% version and documentation can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/sttools/
% Do not use the stfloats baselinefloat ability as IEEE does not allow
% \baselineskip to stretch. Authors submitting work to the IEEE should note
% that IEEE rarely uses double column equations and that authors should try
% to avoid such use. Do not be tempted to use the cuted.sty or midfloat.sty
% packages (also by Sigitas Tolusis) as IEEE does not format its papers in
% such ways.
% Do not attempt to use stfloats with fixltx2e as they are incompatible.
% Instead, use Morten Hogholm'a dblfloatfix which combines the features
% of both fixltx2e and stfloats:
%
% \usepackage{dblfloatfix}
% The latest version can be found at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/dblfloatfix/




%\ifCLASSOPTIONcaptionsoff
%  \usepackage[nomarkers]{endfloat}
% \let\MYoriglatexcaption\caption
% \renewcommand{\caption}[2][\relax]{\MYoriglatexcaption[#2]{#2}}
%\fi
% endfloat.sty was written by James Darrell McCauley, Jeff Goldberg and
% Axel Sommerfeldt. This package may be useful when used in conjunction with
% IEEEtran.cls'  captionsoff option. Some IEEE journals/societies require that
% submissions have lists of figures/tables at the end of the paper and that
% figures/tables without any captions are placed on a page by themselves at
% the end of the document. If needed, the draftcls IEEEtran class option or
% \CLASSINPUTbaselinestretch interface can be used to increase the line
% spacing as well. Be sure and use the nomarkers option of endfloat to
% prevent endfloat from "marking" where the figures would have been placed
% in the text. The two hack lines of code above are a slight modification of
% that suggested by in the endfloat docs (section 8.4.1) to ensure that
% the full captions always appear in the list of figures/tables - even if
% the user used the short optional argument of \caption[]{}.
% IEEE papers do not typically make use of \caption[]'s optional argument,
% so this should not be an issue. A similar trick can be used to disable
% captions of packages such as subfig.sty that lack options to turn off
% the subcaptions:
% For subfig.sty:
% \let\MYorigsubfloat\subfloat
% \renewcommand{\subfloat}[2][\relax]{\MYorigsubfloat[]{#2}}
% However, the above trick will not work if both optional arguments of
% the \subfloat command are used. Furthermore, there needs to be a
% description of each subfigure *somewhere* and endfloat does not add
% subfigure captions to its list of figures. Thus, the best approach is to
% avoid the use of subfigure captions (many IEEE journals avoid them anyway)
% and instead reference/explain all the subfigures within the main caption.
% The latest version of endfloat.sty and its documentation can obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/endfloat/
%
% The IEEEtran \ifCLASSOPTIONcaptionsoff conditional can also be used
% later in the document, say, to conditionally put the References on a
% page by themselves.




% *** PDF, URL AND HYPERLINK PACKAGES ***
%
%\usepackage{url}
% url.sty was written by Donald Arseneau. It provides better support for
% handling and breaking URLs. url.sty is already installed on most LaTeX
% systems. The latest version and documentation can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/url/
% Basically, \url{my_url_here}.




% *** Do not adjust lengths that control margins, column widths, etc. ***
% *** Do not use packages that alter fonts (such as pslatex).         ***
% There should be no need to do such things with IEEEtran.cls V1.6 and later.
% (Unless specifically asked to do so by the journal or conference you plan
% to submit to, of course. )


% correct bad hyphenation here
\hyphenation{op-tical net-works semi-conduc-tor}
%Added by Yilun
\usepackage{amssymb,amsmath}
\usepackage[normalem]{ulem}
\usepackage[color,final]{showkeys}
\usepackage{multirow}
\usepackage{xhfill}
\usepackage{graphics}
\usepackage{graphicx}
\usepackage{subfigure}
\usepackage{epsfig}
\usepackage{url}

\newcommand{\remove}[1]{{\color{Gray}\sout{#1}}}
\newcommand{\comm}[1]{{\color{red}(#1)}}
\newcommand{\revise}[1]{{\color{blue}#1}}
\newcommand{\R}{{\mathbb{R}}}
\def\ni{\noindent}

\begin{document}
%
% paper title
% can use linebreaks \\ within to get better formatting as desired
% Do not put math or special symbols in the title.
\title{Randomized Structural Sparsity  based Support Identification with Applications to  Locating  Activated or Discriminative Brain Areas: A  Multi-center Reproducibility Study}
%
%
% author names and IEEE memberships
% note positions of commas and nonbreaking spaces ( ~ ) LaTeX will not break
% a structure at a ~ so this keeps an author's name from being broken across
% two lines.
% use \thanks{} to gain access to the first footnote area
% a separate \thanks must be used for each paragraph as LaTeX2e's \thanks
% was not built to handle multiple paragraphs
%

\author{Yilun~Wang, %~\IEEEmembership{Member,~IEEE,}
        Sheng~Zhang, %~\IEEEmembership{Fellow,~OSA,}
        Junjie~Zheng,
        Heng~Chen,
        and~Huafu~Chen,  \thanks{This work was supported by the 973 project (Nos. 2015CB856000, 2012CB517901), 863 project (No. 2015AA020505), the Natural Science Foundation of China (Nos. 61125304, 11201054, 91330201), and the Specialized Research Fund for the Doctoral Program of Higher Education of China (No. 20120185110028), and the Fundamental Research Funds for the Central Universities (ZYGX2013Z004, ZYGX2013Z005). Corresponding authors: Yilun Wang and Huafu Chen.}%~\IEEEmembership{Life~Fellow,~IEEE}% <-this % stops a space
\thanks{Yilun~Wang and Sheng Zhang are with the School of Mathematical Sciences  and Center for Information in Biomedicine, University of Electronic Science and Technology of China, Chengdu,
Sichuan, 611731 P. R. China. (email: yilun.wang@rice.edu)
%of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta,
%GA, 30332 USA e-mail: (see http://www.michaelshell.org/contact.html).
}% <-this % stops a space
\thanks{Junjie~Zheng, Heng~Chen and  Huafu~Chen  are with Key laboratory for Neuroinformation of Ministry of Education,
School of Life Science and Technology, and Center for Information in Biomedicine,
University of Electronic Science and Technology of China,
Chengdu, Sichuan, 610054 P. R. China. (email: chenhf@uestc.edu.cn)}% <-this % stops a space
%\thanks{Manuscript received April 19, 2005; revised December 27, 2012.}
}

% note the % following the last \IEEEmembership and also \thanks -
% these prevent an unwanted space from occurring between the last author name
% and the end of the author line. i.e., if you had this:
%
% \author{....lastname \thanks{...} \thanks{...} }
%                     ^------------^------------^----Do not want these spaces!
%
% a space would be appended to the last name and could cause every name on that
% line to be shifted left slightly. This is one of those "LaTeX things". For
% instance, "\textbf{A} \textbf{B}" will typeset as "A B" not "AB". To get
% "AB" then you have to do: "\textbf{A}\textbf{B}"
% \thanks is no different in this regard, so shield the last } of each \thanks
% that ends a line with a % and do not let a space in before the next \thanks.
% Spaces after \IEEEmembership other than the last one are OK (and needed) as
% you are supposed to have spaces between the names. For what it is worth,
% this is a minor point as most people would not even notice if the said evil
% space somehow managed to creep in.



% The paper headers
\markboth{Journal of \LaTeX\ Class Files,~Vol.~11, No.~4, December~2012}%
{Shell \MakeLowercase{\textit{et al.}}: Bare Demo of IEEEtran.cls for Journals}
% The only time the second header will appear is for the odd numbered pages
% after the title page when using the twoside option.
%
% *** Note that you probably will NOT want to include the author's ***
% *** name in the headers of peer review papers.                   ***
% You can use \ifCLASSOPTIONpeerreview for conditional compilation here if
% you desire.




% If you want to put a publisher's ID mark on the page you can do it like
% this:
%\IEEEpubid{0000--0000/00\$00.00~\copyright~2012 IEEE}
% Remember, if you use this you must call \IEEEpubidadjcol in the second
% column for its text to clear the IEEEpubid mark.



% use for special paper notices
%\IEEEspecialpapernotice{(Invited Paper)}




% make the title area
\maketitle

% As a general rule, do not put math, special symbols or citations
% in the abstract or keywords.
\begin{abstract}
In this paper, we focus on how to locate the relevant or discriminative brain regions related with external stimulus or certain mental decease, which is also called support identification, based on the neuroimaging data.
 The main difficulty lies in the extremely high dimensional voxel space and relatively few training samples, easily resulting in an unstable  brain region  discovery (or called feature selection in context of pattern recognition). When the training samples are from different centers and have between-center variations, it will be even harder to obtain a reliable and consistent result. %for many traditional multivarate feature selection methods.
%In order to deal with these difficulties, many efforts have been performed and we will mainly focus on designing stable feature selection algorithms.
Corresponding, we revisit our recently proposed algorithm based on  stability selection and
%
%
%stability selection for discriminative features have recently received great attentions, due to its
%finite sample control for certain error rates of false discoveries and
%transparent principle for choosing a proper amount of regularization.
%However, it fails to consider some prior or estimated structural information of these features. Correspondingly,  we propose a
%algorithm which is able to incorporate the prior
%
structural sparsity.  % is able to  achieve stable  and reliable support identification.  %While its advantages over the common stability selection is the better control of false negatives, yet keeping the control of false positives. %Compared with the common structural sparsity models
%, the prior or estimated structural information is only for a reference and not necessarily very accurate, which is important in practice.
It is applied to the multi-center MRI data analysis  for the first time. A consistent and stable result is achieved across different centers despite the between-center data variation while many other state-of-the-art
methods such as two sample t-test fail. Moreover, we have empirically showed that the performance of this algorithm is  robust and insensitive to several of its key parameters.
In addition, the support identification results on both functional MRI and structural MRI are interpretable and can be  the potential biomarkers.


%In this paper, we focus on the  functional magnetic resonance imaging data of the brain, and  we present a specific implementation based on our proposed
%I.
%
%voxels from functional magnetic resonance imaging data have been gaining great attention recently, because the widely used univariate methods such as the %two-sample t-test fail to find multi-voxel patterns. %, and remain an unsolved problem, because more attention is paid to the prediction accuracy instead.
%In this study, we focus on sparse pattern recognition methods %based identification of the involved brain regions in
%for very high dimensional functional
%MRI (fMRI) data associated with either a cognitive task or disease  diagnosis.
%However, the inferred discriminative voxels  of most existing methods %are instantiated and optimized for the prediction accuracy, but their inferred discriminative voxels
%either have significant false negatives or false positives,
%However, most of them fail to give accurate and stable results in the cases of very high dimensional voxel space. In order to partially overcome this problem, we proposed a  variant of stability selection  based on
%
%``constrained block subsampling" to
%and multi-stage refineness aiming to %it is able to reduce the falsely selected voxels, aiming to
% reliably and stably recover the ground truth of the discriminative voxels  efficiently.
%Numerical experiments
 %We
% demonstrate that our method can achieve  better  accuracy and reproducibility than many alternatives based on the multi-center autism brain imaging data across different modals. %Specifically, compared with the original stability selection, this variant \end{abstract}
% Note that keywords are not normally used for peerreview papers.
\end{abstract}
\begin{IEEEkeywords}
voxel selection; multi-center; structural sparsity; stability selection; constrained block subsampling; fMRI; sMRI; Multimodal;  pattern recognition
\end{IEEEkeywords}






% For peer review papers, you can put extra information on the cover
% page as needed:
% \ifCLASSOPTIONpeerreview
% \begin{center} \bfseries EDICS Category: 3-BBND \end{center}
% \fi
%
% For peerreview papers, this IEEEtran command inserts a page break and
% creates the second title. It will be ignored for other modes.
\IEEEpeerreviewmaketitle



%% main text
\section{Introduction} \label{Sec:intro}
\subsection{Problem Statement}
For computer multimedia analysis, a challenging problem is the semantic gap between the high-level perception and cognition  and the low-level features of the digital contents. Considering that the human brain can eliminate this gap very naturally, deep understanding of the brain responses to multimedia will be playing an important role in designing the computational strategies for multimodal representation, classification and retrieval.
 They share  a fundamental problem with  many other brain related fields such as clinical diagnosis of  mental diseases.  That is how to reliably locate the relevant brain regions corresponding to either external multimedia stimulus or certain mental disease. Via modern brain imaging techniques such as functional magnetic resonance imaging (fMRI), electroencephalography (EEG) and magnetoencephalography (MEG), ones can  image the structure, function or pharmacology of the nervous system in a non-invasive way. Recently, learning from neuroimaging data, as a kind of pattern recognition,  has led to impressive results \cite{Haxby01MVA}, such as  guessing which image or video a subject is looking at
%from his brain activity , as well as medical diagnosis, e.g.,
 or where are the most activated brain regions when one responses to multimedia stimulus in brain disorders. {In many such applications, % besides expecting better prediction scores by designing effective classifiers, we are also trying to  accurately identify the true and reliable discriminative or called predictive features (voxels, in this paper) ,
%it is  quite an  important goal in many application in brain decoding and clinical diagnosis \cite{Li09SparseVoxelSelection,Guyon03Introduction,Blum1997245}. For example,  %because the discovery of biomarkers is typically modeled as selecting the most discriminating and reliable features based on the classification information (e.g. discriminating healthy versus diseased)
the commonly and reproducibly selected brain areas can be considered as the potential biomarkers \cite{Hilario08biomarker,Li09SparseVoxelSelection,Guyon03Introduction,Blum1997245}.  Without reliable, reproducible results, no study in generally can effectively contribute to scientific knowledge \cite{Bennett10reliableFMRI}. }%The reproducibility is quite important for the  Say something about multi-center data analysis.}  Reproducibility and consistency. %In the language of statistics and machine learning, this is often referred to as feature

Discovering the  discriminative brain regions corresponding to  certain mental decease or activated brain regions corresponding to certain
external stimuli, also called support identification,  is a typical  feature selection problem in the context of pattern recognition. %re selection is considered to be  important for better prediction accuracy, though
%It becomes a difficult task when the feature space is of ``high dimensionality" .
%selection. Feature selection has attracted strong research interest in the past several decades. For recent reviews
%of feature selection techniques used in bioinformatics, the reader is referred to [5, 6, 3, 7]
%
%\comm{Multi-center}
%
%
%
However, due to the high dimensionality of feature space and relatively few samples, it is quite a challenging task to accurately and stably estimate discriminative voxels  \cite{buhlmann2011statistics}. Many existing  multi-variable pattern analysis methods often fail to provide a stable feature selection result, i.e., there exists the inherent significant
run to run variability in the decision space generated by the classifiers in terms of  even very small changes to the training set \cite{Anderson10critiqueMVPA}. This instability  is
   partially because classic multivariate feature selection methods aim at selecting a
minimum subset of features to construct a classifier of the best predictive accuracy and  often ignore
``stability" in the algorithm design \cite{Yu08denseFeatureGroups,He2010215,buhlmann2011controlling,Cover65Pattern,Anderson10critiqueMVPA,Raizada2010PFN}. %}
%
%


Therefore, we will focus more on the stability and reliability of the feature selection results, for the sake  of scientific truth. {  Correspondingly, % feature selection methods need to be proposed and their performance needs to verified. Therefore,
%reducing the false positives as well as false negatives during feature selection is an important research topic.
 our research will be based on multi-center neuroimage data analysis because multi-center data allows  to test the
reproducibility and consistency of the feature selection  results. In addition,  a multi-center  study offers several
advantaged over a single center study. For example, it has the potential to increase the number and diversity of subjects enrolled, and to include significant numbers of subjects from
rare clinical subgroups \cite{Mueller10how}. While  inter-center variability brings great challenges for consistent feature selection results, it can also be used to test the performance of different algorithms. %In this paper, we are using disease-related multi-center data, which are commonly used for test the generalization of a feature selection method. Surely, our method can also be extended to support identification in the context of external multimedia stimuli.
Finally, we would like to emphasize that in line with the trend towards Big Data in brain study, major advances have been made in the availability
of shared neuroimaging data. There have existed more than
$8,000$ shared MRI (magnetic resonance imaging) data sets available online \cite{Poldrack14neuroimaging}. The quantitative analysis of these multicenter data sets
 from related magnetic resonance imaging experiments has the potential to significantly accelerate progress in brain
mapping. }


In order to further increase the trustworthiness of the discovered biomarkers, we consider both functional MRI (fMRI, for short) and structural MRI (sMRI, for short). While the results from functional imaging are increasingly being
submitted as evidence into the legal systems of many nations,  %For example, results from a
%commercial company called No Lie MRI (San Diego, CA; http://www.noliemri.com/) were
%introduced into a juvenile sex abuse case in San Diego during the spring of 2009.
%The defense
%was attempting to introduce the fMRI results as scientific justification of their client��s claim of
%innocence. A concerted effort from imaging scientists, including in-person testimony from Marc FMRI RELIABILITY
%Raichle, eventually forced the defense to withdraw the request. While the fMRI results never
%made it into this case, it is clear that fMRI evidence will be increasingly common in the
%courtroom.
the reliability of the evidence from fMRI still needs to further proved \cite{Bennett10reliableFMRI}. Therefore, we also consider the structural MRI data, which is generally considered to be less dynamically changed. The biomarker results from both fMRI and  sMRI will be considered together to give a better and more convincing understanding of the decease, i.e. Autism Spectrum Disorder (ASD) in this paper \cite{Serres14multimodals,Wainer14}.

%\comm{However, feature selection is harder}
%While the prediction or classification accuracy of these designed classifiers are mainly considered in most existing literatures, we are focusing more on the first component, i.e. feature selection, which is aiming to  accurately identify the true and comprehensive  discriminative or called predictive features. This task is also called support identification,  which is even a more important goal in many practical cases including the medical diagnosis  where these selected voxels  can be used as biomarker candidates \cite{Guyon03Introduction}.




%{Therefore, in this paper, we focus on this challenge of data variation between different centers.}
{  To our best knowledge, most of existing algorithms fail to overcome  this difficulty caused by data variation. For example, %As mentioned above, we focus more on the algorithmic point of view to deal with this kind of difficulty arising from sampling variations of MRI data
%viewpoint of statical inference,  the possible unstability of feature selection could be due  to
%gathering  from different sites \cite{He2010215}.
%revisit our prosed ``randomized structural sparsity" based
% A stable feature selection method is expected to be applied  to each of the multi-centers and obtain a consistent feature selection results across different sites.
%
the univariate approaches have been popular for finding the discriminative regions
due to their being directly
testable, easily interpretable, and computationally tractable. However, it is not quite robust to the data perturbation possibly caused by the between-center variation.  In addition,  there exists a strong demand  in various multivariate approaches in
recent years \cite{Haxby14Neural}, because the brain representation
is intrinsically multivariate. However, as mentioned before, many multivariate methods often ignore
``stability" in terms of feature selection. }
% and correspondingly.}
% In other words, non-degenerate classes will always be linearly separable
%\cite{}, and a classifier can then trivially obtain 100\% correct on any training set, without
%it necessarily having accurately captured all  the discriminative  structures of the data \cite{}. In other words, in such cases,
%feature selection is even harder in terms of controlling false positive selections and false negative selections, especially when features are often of strong interrelations \cite{}. Therefore,
%It becomes an even  harder task to accurately and stably perform feature selection in cases of small samples and high dimensional feature space.
 %It is closely related with the we-called ``Laziness" phenomenon of the classifiers in terms of support identification or feature selection, which will be demonstrated in the section of numerical experiments.


  {We aim to achieve a consistent feature selection result across different sites.  Correspondingly,  %In this paper, we introduce the Autism Brain Imaging Data
%Exchange (ABIDE; http://fcon\_1000.projects.nitrc.org/indi/abide/)��
%a grassroots consortium effort dedicated to aggregating and
%sharing previously collected resting-state functional magnetic
%resonance imaging (R-fMRI) data sets from individuals with ASDs
%and age-matched typical controls (TCs).
%
%
%While our methodology can be applied to many other kinds of high dimensional data mining, we would like to take a specific kind of data as an example to show its detailed implementations. In this paper,  we are mainly focusing on the neuroimaing data, where the curse of dimensionality constitutes a major challenge to pattern recognition.  In addition, learning from neuroimaging
%data can lead to many meaningful results, such as brain reading or brain decoding, e.g. guessing which image a subject is looking at
%from his brain activity \cite{Haxby01MVA}, as well as medical diagnosis, e.g., finding out whether the person is a healthy control or a patient.
%While many of the existing pattern recognition works pay more attention to  the prediction or classification accuracy,  we are more focusing on reliable feature selection, i.e. accurately identifying the brain activation regions involved in the certain cognitive processing of an external stimuli or a given disease.
%
%
we will revisit  our recently proposed concept of ``randomized structural sparsity" and  apply it to the multi-center data analysis. We will verify the stability of the resulted algorithm in terms of the generalization and reproducibility  of the selected features via the multi-center data \cite{yu2013}. %At a minimum, reproducibility manifests itself in stability of statistical results relative to ``reasonable" perturbations to data.
In particularly, we will
  %use the multi-center data %we expect the method to be among the most sensitive, i.e. detecting as many true positives as possible. %. \comm{Sensitivity, Receiver operating characteristic, reproducibility }
%\comm{This is related with stability selection}
%
%We will using the multi-center
   %to
   check whether a consistent feature selection result can be obtained across different centers, and  this consistency will support the reliability of the selected features as the potential true biomarker of the associated disease, i.e. ASD here.} %and the effectiveness of our method. \comm{Correspondingly, we will give some potential biomarker to ASD based on multi-center multimodal data.}  %We will verify the stability and reproducibility of the results based on our method based on the multi-center data by considering the consistency of the feature selection results on the different centers. %\comm{Comparison of the predictive ability in terms of sensitivity and specificity of different topological descriptors}
%\comm{Reproducibility: due to the relatively small number of subjects, very high dimensional feature space (consisting of the network edges) and probably high level noise, leading to the over-fitting training data and ����curse of dimensionality���� problem.}

% Each center will be treated individually. Then we will check the consistency of the results of different centers. \comm{While most of previous studies based on multi-center data} mainly consider the , \cite{}, we will focus more on the analysis and verification of our proposed algorithm in \cite{}.

 % mainly consider the function MRI data,   our algorithm might be applied to other kinds of modalities.
%Specifically, from viewpoint of pattern recognition, we consider each fMRI volume (referred
%to as a sample) consisting of $p$ voxels (features) as a vector in a $p$ dimensional feature space and each sample is  labelled or classified reflecting different experimental
%conditions or brain states, and then the obtained relationship between the voxels and their corresponding labels are tested with
%another independent set of samples.


%\comm{Feature Extraction--Maybe removed}The mulivariate methods include supervised approaches such as linear discriminant analysis (LDA), unsupervised ones %such as principal component analysis (PCA), and additional spectral and manifold learning methods. When the class label information is available, supervised
%approaches, such as LDA, are usually more effective than unsupervised ones such as PCA for classification. In addition, the limitation of the feature %extraction methods like PCA is that  the new features created by combining the original features often lack intuitive physical meanings, which hinders result %interpretation.

%verify the reliability of feature selection results of our algorithm by checking their  consistency across different sites.  Clinical and Diagnostic Applications


\subsection{Advantages and Limitations of Sparsity Applied to Nueroimaging}
%\comm{Our focus? Why we only review single site methods.}
{While most of existing efforts to deal with multi-center data is to  merge the data from different sites into a pool  and analyze it as a whole, like a big single center, we  consider each of multiple centers individually \cite{Martino13MP}. } In particular, we consider each of the centers separately and apply the feature selection method to obtain the potential biomarker. Then we  compare the results from different centers to see whether they are consistent to a certain degree. While this seems to be a natural and better way to validate  the reproducibility of the final selected biomarkers than only considering a  single center data, it brings a great challenge to the feature selection method, which is required to be stable due to the possible high data variability across different centers. Now we first review some commonly used feature selection methods in the following parts.


We consider to  identify the discriminative brain voxels from the given labeled training MRI data in context of the supervise learning. % and %do not consider feature extraction and their corresponding methods such as Principle Component Analysis \cite{Jolliffe05PCA}, Linear Discriminant Analysis \cite{Mika99FDA},  Non-negative Matrix Factorization \cite{Lee99NMF}, and  Independent Component Analysis \cite{Hyvarinen:2000:ICA,Daubechies09SparsityOvreICA}.
%As mentioned before,
While the  classification problem is mainly considered, the regression problem can be treated in a similar way.  The linear model is widely  used for feature selection as follows.
  \begin{equation} \label{eq:linearmodel}
 \mathbf{y} = \mathbf{X}\mathbf{w} + \epsilon
 \end{equation}
 where $\mathbf{y} \in \mathbb{R}^{n\times 1}$ is the binary classification  information and $\mathbf{X} \in \mathbb{R}^{n\times p}$ is the given training voxel-wise MRI data and $\mathbf{w} \in \mathbb{R}^{p\times 1}$ is the unknown weights reflecting the
 degree of importance of each voxel.  Identification of discriminative voxels  is based on the values of the weight vector $\mathbf{w}$ and their importance is  usually proportional to the absolute values of the components. %thcan be formulated as detecting the non-zero coefficients of large magnitude, which shows the importance of each voxels in terms of this classification task or certain cognitive processing of an external stimuli.
% In this paper, we %mainly consider the problem of identification of the involved brain regions, i.e. voxel selection, and
%do not consider feature extraction and their corresponding methods such as Principle Component Analysis \cite{Jolliffe05PCA}, Linear Discriminant Analysis %\cite{Mika99FDA},  Non-negative Matrix Factorization \cite{Lee99NMF}, and  Independent Component Analysis \cite{Hyvarinen:2000:ICA,Daubechies09SparsityOvreICA}. %e  our hypothesis on the spatial layout of the neural coding?}

%As mentioned before,
%\comm{$p \gg n$}

Considering that the common challenge in this field is the curse of dimensionality $p \gg n$. %, e.g. the number of
%voxels (tens of thousands, for example) is often much larger than the numbers
%of available (training) samples, for example, hundreds of time points.
%The standard univariate analysis of neuroimaging data could be less sensitive to the size of $p$. However, it could only find some activity regions lack of multivariate patterns.
%when determining the relationship between brain activity and mental states. %In addition, it fails to achieve the prediction goal.
%As for the multivariate methods, among all the supervised statistical learning methods,
%over the last decade, many contributions have proposed based on  statistical learning and more specifically supervised learning methods %\cite{Mitchell04LearningDecodeCognitive}, because %  The recently popular methods are mostly  multivariate, %: the learned classification or prediction function relies on correlations between
%distant brain voxels,
%since
%brain representations are inherently multivariate \cite{Raizada2010PFN}.  While there have existed unsupervised methods such as independent component analysis (ICA) based methods \cite{Dodel20007ICA,Chen2005ICA}, most of the current support identification  methods are based on the multivariate supervised classifiers, i.e.  making use of the known classification information of the training data
%to infer the involved discriminative brain regions. %, because %This presents major difficulties
%to both classificaion and prediction due to the likely overfitting of the training data.
%supervised machine learning techniques are of % being
%increasingly used in the analysis of brain imaging data
%their inherent ability to deal with multi-variate
%data, higher sensibility and possibility of incorporating
%specific prior-information.
%
%In this paper,
%we are focusing on the sparsity based voxel selection methods in this paper, because %Least
%Squares or Ridge Regression \cite{}, are incapable of effectively discriminating
%which areas of the brain mostly contribute to the model��s predictions, in the sense that all voxels
%are expected to contribute to generate a predictive function.
 sparsity makes much sense that  the most discriminative activated voxels are only
a small portion of the total brain voxels \cite{Yamashita08SparseEstimation}. However, %, even though the majority of these voxels  are involved in the cognitive task or certain disease.  %\comm{},% activated voxels are usually very small in number and are
%embedded in a mass of inactive voxels.
%
%
%We first review some commonly used patter recognition algorithms which can be used for voxel selection. %In such cases,  the label information for each sample %is available and most of corresponding methods are supervised ones.
% On the contrary, the mulivariate methods including supervised approaches such as Support Vector Machine and Logistic Regression can achieve simultaneous identification %of discriminative regions and effective prediction.  The labeled data of samples are used as the training data and each voxel will be corresponding to a calculated %coefficient.
 %Identification of involved voxels can be formulated as detecting the non-zero coefficients of large magnitude, which shows the importance of each voxels in terms of %this classification task or certain cognitive processing of an external stimuli.
%\comm{The typical procedure--training and prediction} \comm{Through weights of each voxel, %we get the biomarkers.} Most of the existing multivariate pattern recognition algorithms focus more on the prediction or classification accuracy than the %identification of the meaningful biomarkers.
%
%Due to the fact that sparsity makes us reconstruct high dimensional data with few samples,
%
%For better support identification performance, the estimators should incorporate this prior knowledge. It is natural to promote prediction
%functions relying on only a few brain regions, for instance using sparsity
%inducing regularization methods  for example, the well known $\ell_1$ norm regularization.   Due to the fact that sparsity makes us reconstruct high dimensional data %with few samples,  here we are focusing on sparsity based feature selection methods, which is also more relevant than independence as identification criterion than ICA %\cite{Daubechies09SparsityOvreICA}.  After all, only a small fraction of the brain is involved in the task and this is where the sparsity comes from.
sparsity alone is not
sufficient for making reasonable and stable inferences, in cases of very high voxel space and small number of training samples. In such cases, the plain sparse learning models such as the $\ell_1$-norm regularized model also called LASSO \cite{tibshirani1996regression},  often provide
overly sparse and hard to interpret solutions %where % where
%the non-zero coefficients or
%the  selected voxels are often assigned to be individually scattered
% disparateregions across the brain  in our cases
\cite{Rasmussen2012SparsityFMRI}. % with either high false positive rate or false negative rate. %due to the very small number of samples.
%As for the fMRI,  the data of neighboring voxels might be of strong correlation and the plain $\ell_1$ norm based models  tend to
%pick only one or a few of them and shrinks the rest to $0$. %This may not be
%a desirable feature. For example, in microarray analysis, expression levels of
%genes that share one common biological pathway are usually highly correlated, and these genes may all contribute to the biological process, but lasso
%In addition, things can be even worse, due to the the strong local correlation between columns of $\mathbf{X}$,  the model (\ref{eq:L1}) itself could  %be find the involved brain regions effectively in an interpretable way, because
 %likely makes  the plain $\ell_1$ model
%  select arbitrarily one of them and not the others,  resulting into  a high rate of false negatives in the
%support identification.  In such cases,
Thus we need to extend the plain sparse learning model to incorporate some important structural features of brain imaging data in order to achieve a more stable, reliable and interpretable support identification results.




%the grouping or correlation information of neighboring or contiguous voxels of fMRI data.
%Thus we might extend the plain $\ell_1$ model to a more flexible structure sparsity induced regularization to incorporate the grouping or correlation information of neighboring or contiguous voxels of fMRI data. %\comm{In other words, overcome the redundancy in terms of classifications}



%because  some important features such as brain segregation and integration, not being fully exploited.
%Functional MRI (fMRI) gives images of brain activity
%via Blood Oxygen-Level Dependent (BOLD) signal changes, and  neuroimaging based on fMRI  is a spatially-smoothed representation
%of the underlying neural signals.
%
%\comm{Specifically,} the neuroimaging data we are focusing on in this paper is the Functional MRI (fMRI), which is based on  Blood Oxygen-Level Dependent (BOLD) signal contrast to gives images of brain activity, though our algorithm might be applied to other kinds of modalities. It is the latest in a line of techniques to study human brain function
%that use brain perfusion as a proxy for brain function (\cite{Pekar06fMRI}),  instead of making any direct measurements of brain function. An fMRI scanner measures
%the blood-oxygenation-level-dependent (BOLD) signal at all
%points in a $3$-D grid, and each grid of the 3-D
%image is known as a voxel. A typical fMRI dataset is composed
%of time series (BOLD signals) of thousands of voxels.



%BOLD fMRI relies on detecting perfusion changes in a brain
%region following a neuronal activity and  neuroimaging based on fMRI  is a spatially-smoothed representation
%of the underlying neural signals.
%Without making use of these features, %Due to the number of the involved regions being much smaller than that of the involved voxels,
%this characteristic  makes
%it is possible to
%select more variables than the number of observations,
%which would be impossible with a standard $\ell_1$-type
%estimator.  \comm{The feature of fMRI data}%\comm{Localization corresponding to strong Spacial Correlation. }
%
%\comm{Learning with incomplete data}
% While sparsity has traditionally been connected
%with interpretability, in the sense that sparser models
%are easier to interpret,
%However, the popular sparsity regularization based on the plain $\ell_1$ minimization (\ref{eq:L1}) has encountered inherent difficulties when applied to this support %identification problem, due to very high dimensional voxel space as well as some important features such as brain segregation and integration, not being fully %exploited. Specifically, %while the simple $\ell_1$ minimization based sparse methods, like the LASSO \cite{tibshirani1996regression} or (\ref{eq:L1}) is able to estimate solutions for which
%only few voxels are deemed relevant, %therefore aiding
%interpretation in some cases,
%the plain sparse learning models often provide
%overly sparse and hard to interpret solutions where % where
%the non-zero coefficients or
%the  selected voxels are often assigned to be individually scattered
%% disparateregions across the brain  in our cases
%\cite{Rasmussen2012SparsityFMRI}. %due to the very small number of samples.
%As for the fMRI,  the data of neighboring voxels might be of strong correlation and the plain $\ell_1$ norm based models  tend to
%pick only one or a few of them and shrinks the rest to $0$. %This may not be
%a desirable feature. For example, in microarray analysis, expression levels of
%genes that share one common biological pathway are usually highly correlated, and these genes may all contribute to the biological process, but lasso
%In addition, things can be even worse, due to the the strong local correlation between columns of $\mathbf{X}$,  the model (\ref{eq:L1}) itself could  %be find the involved brain regions effectively in an interpretable way, because
 %likely makes  the plain $\ell_1$ model
%  select arbitrarily one of them and not the others,  resulting into  a high rate of false negatives in the
%support identification.  In such cases,
%Thus we might extend the plain sparse learning model to incorporate the grouping or correlation information of neighboring or contiguous voxels of fMRI data.
%Thus we might extend the plain $\ell_1$ model to a more flexible structure sparsity induced regularization to incorporate the grouping or correlation information of neighboring or contiguous voxels of fMRI data. %\comm{In other words, overcome the redundancy in terms of classifications}
%\comm{Structures=correlations. randomized for stability. Need to be revised in more unified way.}

%different things. The previous one is mainly about the curse of dimensionality. %This point will also be  demonstrated in our numerical experiments.

 % In addition, randomlization is also very helpful in cases of high dimensional voxel space for stable and reliable results. %and based resampling or subsamping. We will give a brief review of them in the following parts.} %like to combine both in a more efficient way }



%Correspondingly,  we aim to propose an efficient algorithmic framework of sparsity based voxel selection models %to overcome the above difficulty and achieve best performance
%suitable for brain mapping based on the nueroimage  data, with the purpose to  increase sensitivity under a tight control of   false positives, by making use of the above grouping and localisation assumptions.
%


%\comm{not accurate.}  %balance  between computational efficiency and statistical inference quality is also expected to achieve (\cite{Chandrasekaran26032013}). % and ones mostly focus on the last one.  %Here ``false positive" means that we wrongly take some nondiscriminatory voxels as the true discriminative ones. ``false negative" means that we wrongly missed some true discriminative voxels.  \comm{Sensitivity } %both effectiveness and efficiency.
\subsection{Existing  Extensions of the Plain Sparse Model}
% for it is unlikely that only such a small number of genes are involved in the development of a complex disease. A method
%that is able to identify more than n variables should be more desirable for
%such problems.}

Functional segregation of local brain areas that differ in their anatomy and physiology contrasts sharply with their global integration during perception and behavior \cite{Tononi94segregation,Xiang:2009}, and   %MRI responses consisting of voxels whose three dimensional spatial arrangement respect the anatomy of the brain and the discriminative voxels are
%thus expected to have a specific localized spatial organization, which is important for the subsequent identification task performed by neuroscientists .  That is to say,
sets of discriminative or activated voxels between different brain states are expected to form distributed and localized  areas \cite{Chu02featureselection,Tononi94segregation,LimbrickOldfield20121230}.  %We can obtain this kind of structural information either from the anatomical  homology or data driven clustering \cite{HBM:HBM21333,Pekar06fMRI}. % though \comm{not accurate.} %  though the isolated ones are not necessarily false positives, but hard to tell. %un %\comm{Why correlated?}
%
 %the underlying structure needs to be estimated by, for example, some clustering algorithm, and might %be biased.
%
%as the natural supplement of sparsity in order to
%obtain interpretable inferences from almost any sparsity based supervised
%learning methods.
%
%
%
%\begin{itemize}
%\item the role of clustering?  \comm{Unsupervised Dimensional reduction and reduce the correlation.}
%\end{itemize}
 %We take the functional MRI data, which is a kind of temporal-spacial data as the example, though our methods might be extended to other kinds of brain images.
%
%
%\section{}
%\subsection{Review of SPM\_GLM}
%Subsection text here.
%
%
%
%\comm{Why we do not consider the prediction accuracy?}
%\begin{itemize}
%\item
%\comm{Structure of the resulting maps of weights}
%\item How to determine whether our method is the best or not?
%\end{itemize}
%
 %he advantage of our algorithm over the traditional structured sparsity   is better efficiency}
%
%We review some existing  sparsity based multi-voxel support identification works, which   go beyond the plain $\ell_1$ models,  and briefly summarize their advantages %and limitations.
%In such cases, there are two ways to obtain the involved regions, which usually consists of clustered voxels. One is to make use of clustering and the other is turning %to randomlization such as stability selection. %In this paper, We would like to combine both in a more efficient way }
% which is only partially solving the above problems. }
%\begin{itemize}
%\item
%\item Challenge: the true model does not belong to $\ell_1$ regularization path Approach
%\item Stabilize the sparse estimator with randomization [Meinshausen and
%B��hlman, 2009]
%
%\item We are not only interested in the  prediction or classification, but also on brain
%mapping, that is identifying the brain regions involved
%in the cognitive processing of an external stimuli or a given disease.  The best predictive model
%may have selected by chance non-informative regions, and neglected relevant
%regions that provide duplicate information. So we need to
%address the support identification problem \cite{Mitchell04LearningDecodeCognitive}.
%\item Another challenge of high-dimensional sparsity recovery
%from noisy observations is the difficulty to control
%false positives,
%\end{itemize}
%\comm{Functional MRI data are a spatially-smoothed representation of the underlying neural signals. Consequently, the activations are spatially correlated}
Therefore,  two common hypothesis have been made for MRI data analysis. One is the sparsity: the discriminative or activated voxels implied in the classification task or an external stimulus are only a small portion of the total; the other is
 compact structure: these voxels are often grouped into several distributed  clusters because those within a cluster have similar behaviors. % and correspondingly strongly  correlated. %for fMRI data, the very high-dimensional voxel space makes simple sparsity model based on the pure $\ell_1$ penalty is far away from enough, and often results in both false negatives and false positives \cite{Li09VoxelSelection}.
%The brain is of regional segregation and  correspondingly the fMRI data are often highly correlated spatially and the discriminative voxels often consist of several distributed clusters of internal multivariate voxel interactions.  While the above fact might brings computational difficulties due to the ill-conditioning  of $\mathbf{X}$,  it also provides ones with the possibility to recover the true sets of discriminative voxels, whose number can be much larger than  $n$.
%Thus making use of the above two hypothesis is very important and we will review some state of art existing works along this direction,  %}  %discriminative .  %\comm{Correlation are both bad and good things (clustering for dimensional reduction )}.
 %considering that t
 Notice that the plain sparse learning model based on  the $\ell_1$ regularization is only making use of the first sparsity hypothesis and discard the structure related assumption. %and there have existed some efforts beyond it, as follows.

%As mentioned above, if there is a group of highly correlated variables, then the pure $\ell_1$ regularization model tends to select one variable from a group %and ignore the others.
Elastic Net \cite{Zou05ElsticNet} tries to make use of the voxel correlation by adding an $\ell_2$ regularization or called Tikhonov regularization to the classical $\ell_1$ penalty. The strict convexity due to the help of the added $\ell_2$ term helps to select all together a group of redundant features or none of them in what is known as the grouping effect \cite{RyaliCSM12L1L2Network,Ryali10ElsticNetLogistic}. %term helps % for establishing a network, and  % While for highly correlated features, LASSO tends to select a single feature to represent the %entire group, Elastic net,
%, because  However, the task of balancing the $\ell_1$ and $\ell_2$ terms is not trivial.
%Recently,  other penalties have been added to consider the correlated features besides the Tikhonov regularization
 {Recently, Total-Variation (TV) penalization was added to both $\ell_1$ penalization for voxel selection  \cite{Gramfort:2013:IPR:2552484.2552506}. The TV penalization is used to make use of the assumption that the activations are spatially correlated and the weights of the voxels are close to piece-wise constant. } %Elastic Net has been used for feature selection or pattern classification based on the fMRI data \cite{RyaliCSM12L1L2Network,Ryali10ElsticNetLogistic} and has achieved promising results. %While the Elastic net might achieve improved  sensitivity  than the plain $\ell_1$ model, it does not need or fails to make use any explicit prior knowledge of grouping.  %explicitly made use of the localized spacial contiguousness of discriminative voxels.  \comm{For that reason one may not
%need any prior knowledge of grouping.} %  not be able to identify all the involved most discriminative brain regions. %, though each of those  obtained  discriminative regions is spatially contiguous, which fits
%well with our prior knowledge about the brain��s specialized
%regions mentioned above, as showed in the numerical experiments of this paper.  %brain regions, it might not    %\comm{Elastic net is a related technique. Use elastic net when you have several highly correlated variables. lasso provides elastic net regularization when you set the Alpha name-value pair to a number strictly between 0 and 1.}



 %In this paper, we make our hypothesis on the spatial layout of the neural coding: sparsity, i.e. few relevant voxels/regions implied in the cognitive task, %and compact structure, i.e. relevant voxels grouped into connected clusters.
% Correspondingly we mainly focusing on  such block structured sparsity, i.e. the brain activated or involved voxels are often clustered into several possiblly distributed regions. Notice that we often do not have the information about the grouping beforehand.
%might be time-consuming, or even prohibitive % for real large scale brain decoding problems, due to at least the following $3$ reason.
 %\comm{the advantage of our algorithm over the traditional structured sparsity   is better efficiency}


Furthermore, structured sparsity models are proposed  to  more explicitly  make use of the segregation and integration of the brain  \cite{Bach12StructuredSparsity}, in order to extend the well-known plain $\ell_1$ models. They enforce more structured
constraints on the solution, %and  these new structured sparsity
%models promise an even greater ease of interpretation
%of the identified discriminative patterns, because
such as the discriminative voxels are
grouped together in possibly few clusters \cite{Baldassarre12StructuredSparsityFMRI,Michael11TV,Ye12groupNonconvex,Liu10NIPSGroupeTree,Yuan13Overlapping,Jacob09GLO,Ng11classifier,Ye09L21}. %,  That is to say, existing relationships and structures between
%the voxels (e.g., spatial, hierarchical or related to the physics of the problem at hand)
%are not disregarded, i.e. such specific spatial constraints are incorporated in order for better interpretability.
%Among them, one of the most popular structured sparsity models is the group sparsity regularization \cite{}, where the (possibly overlapping) groups have often  been known as a prior information. % because  in many applications, some intrinsic structure information is available on either input or output sides.
%    For example, multi-task feature learning via efficient $\ell_{2,1}$-norm minimization  has been proposed in \cite{Ye09L21}, and the group sparsity based classifiers has been  applied to fMRI classification in (\cite{Ng11classifier}).
In many cases, the parcellation information is not available beforehand, and therefore, ones can either use the anatomical regions as an approx \cite{Batmanghelich12Groupsparisty}, or  use the data driven methods to obtain the grouping information %such as  the hierarchical agglomerative clustering (Ward hierarchical clustering, for example) and a top-down step to prune the generated tree of hierarchical clusters for obtaining the grouping information
\cite{Michel2012SupervisedClustering,jenatton12multiscale}.
Notice that the kind of parcellation based on certain definition of homogeneousness  may not exactly match the ground truth of the discriminative voxels. However, most structural sparsity models tend to select either all voxels of a group or none of them, and fail to perform further refinements by selecting part of voxels of a group. %andpartially due to not making use of the labelling information $\mathbf{y}$ of the training data \cite{LimbrickOldfield20121230}. %For example, given a parcel, it is possible that only a subregion of it is the discriminative region, not all of it .

% ROI analyses rely strongly on a-priori predictions, so may miss unexpected results. Furthermore, the inclusion of non-active voxels in the ROI average can remove genuine effects, and increase the likelihood of accepting false negatives. }\ % In \cite{} also shared the similar ideas but focus more on the hierarchical sparsity-inducing norms. % In \cite{Gramfort12RandomizedSparsityClustering} adopts the spatially-constrained agglomerative clustering (Ward hierarchical clustering, here) for reduced feature space and sparse logistic regression (based on $\ell_1$ norm) with aiming to achieves both improved prediction and identification for the binary classification problem.
%

%
%However, the structured sparsity can not thoroughly deal with very high dimensional feature space when the number of samples is very small. The reason is that
%So the results  might be biased by the use of certain single grouping setting.

Another important kind of feature selection methods for high dimensional data analysis is stability selection \cite{Meinshausen10StabilitySelection,Shah13StabilitySelection,Hastie09Ensemle}, which %as  a special ``ensemble
%learning" procedures following \cite{},
is  based on resamplings (bootstrapping would behave similarly) and  %In particular, randomization of stability selection  in terms of rows of $\mathbf{X}$, i.e. samples,  randomly takes only a fraction of the training samples and  the columns of $\mathbf{X}$ are reweigthed by a random weight uniformaly sampled in $[\alpha,1]$, where $\alpha$ is a positive number.
%Here, stability selection, in terms of the model (\ref{eq:L1}), is repeatedly running it, except that $\mathbf{X}$ is resampled and reweighted.
%However,
aims to alleviate the disadvantage of the plain $\ell_1$ model, which either selected by chance non-informative regions, or neglected relevant
regions that provide duplicate or redundant classification information \cite{Mitchell04LearningDecodeCognitive,Li12SparseLocalization,Anderson10critiqueMVPA,Poldrack06Criti,Ye12SparseLearningSS,Cao2014Sparse,Ryali20123852}.
%partially because of the
%worrying instability and potential deceptiveness of the most
%informative voxel sets when information is non-local or
%distributed \cite{Anderson10critiqueMVPA,Poldrack06Criti}.
One main advantages of stability selection is the finite sample control of false positives. It also makes the  choice of sparsity penalty parameter  much less critical or sensitive \cite{buhlmann2011statistics}. However, it fails to explicitly make use of
the assumption that these targeted voxels are often spatially contiguous and form into distributed brain regions.
{Correspondingly,  we proposed so called ``randomized structural sparsity" based feature selction algorithm in \cite{Wang14RSS}, by incorporating the idea of structural sparsity into stability selection.  %While it has demonstrate its advantages over many other state-of-the-art algorithms in terms of both better control of false positive rate  and false negative rate on a single site, its performance on multi-center data has not been investigated.
}

% The stability selection  has been applied to the pattern recognition based on brain MRI data and achieved better results than the plain $\ell_1$ models  \cite{}.
%
%
%\comm{%Why Stability Selection? In addition, like other existing  multi-variate classifiers, it %\comm{However,} the existing research has showed that the  multi-variate classifiers
%may
%The latter one is common and more serious for medical diagnosis, and we name it as the the "laziness" of classifiers in terms of support identification,
%While %Given the high-dimensionality of neuroimaging,
%and especially fMRI, data and the few number of samples,
%most of current multi-variate classifiers including supervised approaches such as Support Vector Machine and Logistic Regression
%though  a satisfying prediction accuracy might be still achieved }. % they are not surely to be capable of discriminating all the areas of the brain mostly contributing
%to the model's predictions ,
%.
%Specifically, we demonstrate the cases where %need to overcome a common difficulty of the existing supervised classifiers based support identification methods, i.e.
%only partial discriminative regions can be recovered which we name as , though  a satisfying prediction accuracy might be achieved.
%This ``laziness" phenomena is related with the so-called peaking effect described in \cite{Hughes68} and we will demonstrate it through synthetic data.}  It is also just related with false positives. % But we will consider
%Therefore,  it is expect to figure out  a new effective way for better support identification. %, which means finding out the comprehensive sets of related brain regions, not just parts of them.}
%
%
%\comm{The relationship between the above two.}
%
%While the stability selection is compatible with  the possibly strongly correlated discriminative voxels, it has not made use of the assumption that
%In order to make use of

%However, it is usually computationally prohibitive in order to obtain these locally spatially contiguous discriminative regions, because too many times of resamplings are required due to the very high dimensional voxel space.  Without enough resamplings, stability selection fail to produce  continuous discriminative regions. %\comm{It has some extra drawbacks.}  %Why? In more details}
%
  %Correspondingly, an \comm{exploratory} method is required to accurately identify the truly grouping information of the discriminative voxels. %While this procedure While a high prediction accuracy might be still achieved, the identified brain map might be not correct, or comprehensive or stable though they are indeed cluster-structured and easily readable compared to the overly sparse patterns by the single $\ell_1$ norm minimization models.  In addition, though the incorporation of grouping information, the result might be still not stable due to the too high dimension of the voxel space.  In such cases, the randomization techniques have been adopted  for stable and consistent support recovery \cite{Meinshausen10StabilitySelection}.  %However, randomization itself can only partially solve this kind problem.
%We needs the randomness together with the clusterings, especially for high dimensional correlated feature spaces.
%In addition, the resampling technique is not used in this paper maybe because the
%  computational cost of solving a hierarchical structured sparse optimization problem is big and the resampling is not affordable.
%
%\comm{ How to deal with the grouping:  When and why the "groups" or "structures" are known?}
%
%
%In order to overcome this non-contiguous problem,
%ones  proposed to use the idea of common stability selection together with clustering \cite{Gramfort12RandomizedSparsityClustering,Varoquaux12Clustering}. However, their implementations required running a clustering during each resampling  and %correspondingly the
%computational cost is still almost prohibitive for very large scale problems, due to the high computational cost of the general clustering algorithms. More important,
%
%their adopted ``rescaling" for the columns of $\mathbf{X}$ fails to consider the second hypothesis, i.e. the random reweighting  during their %we do not think it is an efficient way to explore all the possible grouping sets, because
%   implemented stability selection %their rescaling scheme %as for the mining the involved clustered voxels in terms of the labeling information, because the prior reweighting or scaling before clustering
% is voxel-wise and fails to consider the spacial contiguousness of the clustered discriminative voxels. Moreover, their evaluation of different methods are more on the prediction accuracy or classification accuracy, instead of the accuracy of the detected support.
 % Therefore, we need to turn to new resampling implementation for better efficiency \cite{Huang2011LSS}.
%.  \comm{Data driven clustering is based on the original data, not the randomized data. } %Specifically,
%in randomly perturbing the design matrix $\mathbf{X}$ by taking
%only a fraction of the training samples and randomly
%rescaling each variable, in our case each column of $\mathbf{X}$. Then  the clustering algorithm is run and this way
%clustering result is different for every
%randomization, as was done in the Random
%Forests algorithm \cite{Breiman01RandomForests}. And the resulted models such as (\ref{eq:structuralL1}) or its variant are also different.
%By repeating the later procedure and then counting how often each voxel is estimated to be involved, we can get a
%score for each voxel.
% trouble of bias of clustering.  While its % gave a comprehensive presentation of the similar idea in terms of sparse linear regression applied to both regression and classification.  In all the above the resampling techique is adopted due to the instability of the estimated
%support of w, particularly when the columns of X are very correlated as it is the
%case with fMRI data.
%Because of the randomness of the cluster shapes
%in the different iterations, the selection frequency score
%can outline the true discriminative regions, which are even   of a smaller size than the cluster
%size.


%\comm{Since the clustering could be erroneous and highly depends on the number of clusters, we need to take these into consideration when designing an algorithm. One way is to introduce variability through perturbing the clustering. BUT we do think we can do it in a better way via Block Subsamplings.}



%In order to obtain the clustered or locally contiguous results, % overcome the problem of randomoized $\ell_1$ in terms of too many resamplings due to the failure to make use of the prior information about clustering property of the significant discriminative voxels,
% \cite{Varoquaux12Clustering} propose to incorporate the idea of stability selection with clustering.
% However, their method is still quite time-consuming, considering that for each resampling, the adopted clustering algorithm such as the hierarchical clustering algorithm \cite{Murtagh85} needs to be run once, which is usually taking a significant time for high dimensional voxel space. %Thus, and many times of this kind of randomized clustering needs to be run  during the procedure of stability selection.
% In addition, the random reweighting  during their %we do not think it is an efficient way to explore all the possible grouping sets, because
 %  implemented stability selection %their rescaling scheme %as for the mining the involved clustered voxels in terms of the labeling information, because the prior reweighting or scaling before clustering
% is voxel-wise and fails to consider the clustering feature of the neighboring discriminative voxels.
%  Therefore, we need to turn to new resampling implementation for better efficiency \cite{Huang2011LSS}.

%it indeed opens a door to solve the widely known unreliable support identification problems of multivariate classifiers beyond the plain stability selection ,
%and we will give a brief introduction to it later. %because we need to run the hierarchical clustering once in each step of the resamping procedure. %} %because it tends to produce an ��average�� tree which balances with the greedy hierarchical construction of a single tree.  \comm{Solved a reduced %dimensional L1. }
 %The two components are how to cluster and how to define (train and test) the classification model.

%\paragraph{}
%\comm{Constrained Block Subsamplings can bring even better false positive rate.}



%\comm{Network}\cite{Bellec06ClusteringNetwork}, clustering of voxels of homogeneous temporal activity.

%\comm{The difference with other alternative works} Not as before, the sparse linear (logistic)  regression is the only main tool %\cite{Gramfort12RandomizedSparsityClustering,jenatton12multiscale,Varoquaux12Clustering}, our contribution is to incorporate SVM (or LDA \cite{Sato09EvaluateSVMMLDA}) %and GLM into our framework where the sparse  linear (logistic)  regression is only used for supervised feature selection. SVM is for both support identification and %classification. GLM is to
%extend the subset of the discriminative regions.
%\comm{After clustering, Other multivariate methods such as linear discriminant analysis (LDA)
% or SVM   are also capable of
%serving as multivariate ranking methods as . }


 %in more details,  particularly via the adoption of  the ``constrained  block subsampling" technique, instead of single voxel-wise subsampling,  for the selection of discriminative voxels for the fMRI data. %Compared the work in \cite{Varoquaux12Clustering},  to replace its adopted randomized rescaled clusterings,
 %The ``block subsampling" aims to make use of the spatially contiguousness of  discriminative voxels. % in a computationally efficient way. \comm{Empirically it has better control of false positive rate or better balance between false positive rate and false negative rate.} \comm{ Intuitively, correlated voxels are often neighbors, resulting into localized regions. }  %Compared with the classical randomization based stability selection without block %subsamping which mainly aims to reduce the false positive rate, our new method can also reduce the false negative rate.
 %This idea can be easily applied to many existing multivariate classifiers such as sparse SVM, sparse logistic regression, Elastic Net and so on. For the same %inference quality, it is much faster than the common stability selection.
 %The ``constrained" means that we propose to make use of the prior parcelling information of brain or data driven parcelling information to separate or combined the obtained blocked subsamples \cite{}.

\subsection{Our focus and  contributions}
%can help verify the reliability of the selected biomarkers and help remove the false postives}, and a stable feature selection method is required.

%\comm{
  %\comm{Usually computationally prohibitive, Why? In general}
%
%While the high classification or prediction accuracy might be achieved, most of existing algorithms based on multivariate classifiers fail to find out the true %discriminative brain regions, partially because they  are focusing too much on the prediction accuracy, not the identification of related brain regions %\cite{Rasmussen2012SparsityFMRI}. These two goals are not necessarily the same in many cases. %Due to the very small number of training and test samples,
%The smaller number of selected features, i.e. discriminative voxels might still achieve even better prediction accuracy, or at least not worse prediction or %classification accuracy. In other words, the current multivariate classifiers is "active" for good prediction accuracy, but "lazy" in recovering discriminative %voxels, due to the existence of redundant feature sets that can lead to different models of equal predictive power.

% for the very limited number of testing samples.
%The other reason  is due to the very high dimension of the voxel space and the strong correlation between some neighbouring voxels, randomization is required for stable %and reliable and comprehensive support identification. However, most of existing methods fail to make use of it, or not use it in a computationally efficient way. % adoption of resampling techniques for stability, accuracy, and completeness of the identified brain regions. The resampling is not only for numerical stability, but also for the gradually outlining the relatively accurate boundaries of the involved brain regions by running different trials of subproblems. Without making proper use of the specific features of the set of the discriminative voxels, too many times of resamplings are required, especially for the extreme high dimensional problem and it might be computationally prohibitive.
%The existing work performs the "guessing" by the clustering, which is quite computationally costly in high dimensional voxel spaces \cite{Varoquaux12Clustering}. % in order

%\item The quality of spatial patterns extracted from models cannot be assessed purely by focusing on prediction accuracy .
%\end{itemize}

%\comm{in more details. the advantage of our algorithm is better efficiency}
%\comm{The important role of randomness?}

%\comm{The analysis approach taken for multi-center imaging data will depend on the size of the task activation and the
%amount of residual site effect observed in the data. If the between-site reliability is known to be high (say 0.75 or
%higher), it may be acceptable ignore site in the analysis, simply merging the data into a single pool and performing
%the analysis.}

%The key point of this paper is to  achieve both computational efficiency and better inferential quality in terms of both smaller false positive rate and false negative rate, by making proper use of the above two hypothesis, especially that correlated discriminative voxels are often neighbored and locally clustered.  In addition,  stability via subsampling or
%bootstrapping which has been applied in the context of
%prediction \cite{Breiman96Bagging} as well as for the variable selection \cite{buhlmann2011statistics}, will be incorporated for further improving the stability of the feature selection result.
%We aim to achieve a better balance between controlling for false positives and controlling for false negatives. %ational efficiency on the common computing platform. %\comm{Correlated voxels are also neighbored. }
% Correspondingly,  we are aiming to propose a computationally efficient  algorithmic framework  based on the tailored generalization of the classical stability selection. %in a computationally efficient way. % make use of the idea of stability selection \cite{} and first propose make use of blocked resampling on the voxels.}
%\comm{While the stability selection is using the resamping techniques, which means repeating solving a sequence of sparse optimization problems, but the size %of each such problem is small.}
%
%\comm{Mathematically} Our further development is the balance of the effectiveness and efficiency.
 % We are focusing the methodology based on the resampling, clustering and sparse optimization.
%
% Correspondingly,  our main methodological contribution is to propose a general framework for efficient brain mapping, which is usually large scale very high %dimensional problems. Specifically,  %but adopt some simple but faster clustering algorithms such as K-means algorithms for computational efficiency for large scale problems. At the first, we can even just using simple partition techniques at the beginning of our algorithms.
{In this paper, we consider to achieve a consistent and reproducible  feature selection result.  %separately, it
 We revisit the ``randomized structural sparsity" based feature selection method we recently proposed in \cite{Wang14RSS}, where the finite sample control of false positives and false negatives are considered on a single-center study. Here we demonstrate its reproducibility and stability of the selected biomarkers verified through the multi-center data analysis. % compared with  the state-of-the-art methods such as two sample t-test. SVM, logistic regression  \cite{Carroll2009112,Charalampous12sparseDeepLearning,Michael11TV}. % particularly via the adoption of  the ``constrained  block subsampling" technique, instead of single voxel-wise subsampling,  for the selection of discriminative voxels for the fMRI data. %Compared the work in \cite{Varoquaux12Clustering},  to replace its adopted randomized rescaled clusterings,
% The ``block subsampling" aims to make use of the spatially contiguousness of  discriminative voxels. % in a computationally efficient way. \comm{Empirically it has better control of false positive rate or better balance between false positive rate and false negative rate.} \comm{ Intuitively, correlated voxels are often neighbors, resulting into localized regions. }  %Compared with the classical randomization based stability selection without block %subsamping which mainly aims to reduce the false positive rate, our new method can also reduce the false negative rate.
 %This idea can be easily applied to many existing multivariate classifiers such as sparse SVM, sparse logistic regression, Elastic Net and so on. For the same %inference quality, it is much faster than the common stability selection.
 %The ``constrained" means that we propose to make use of the prior parcelling information of brain or data driven parcelling information to separate or combined the obtained blocked subsamples.  %to guide or constrain the block subsampling, to further accelerate the convergence of our stability selection, i.e. reducing the number of resemblings, as well as improving the empirical selection accuracy. \comm{Together with data driven analysis}% of empirical  . %, namely constrained block subsmapling. Compared with the %it is a computationally efficient way to explore most likely groupings among all possible sets of grouping to   quickly find out  all the true discriminative brain regions, which are clusters of voxels.
 %Thirdly, \comm{where ��?>?0 is a given threshold value, i.e., a feature is finally selected if, for at least one value of ��, the fraction of bootstra p experiments for %which the feature is selected exceeds the threshold ��.}
 % Thirdly, for support identification, we do not have to use the group-sparsity induced norm as most existing works have done. Instead we only use the average of each cluster of  voxels and solve a reduced size plain $\ell_1$ model. The reason is that we only focus on the support information in this paper, rather than the accurate weighting value for each voxels. In addition, we can use the existing $\ell_1$ solvers and do not have to develop new group sparsity solvers. Also, due to the reduced dimension, we expect to achieve a much higher  computational efficiency and more stable estimation. %, because the averaging is playing an important role in dimensional reduction like PCA.  % much cheaper way.
%Third unlike the existing preliminary works, we are trying to reduce the number of clustering, which is time-consuming and occupy most of the total running time.
%Third,  we are is to make positive use of the correlations
%in the design matrix, to  improve the behavior of sparse methods, rather than passively taking the correlations as negative. Our
%method can thus select more voxels than the number of observations, which
%would be impossible with standard SLR and difficult with only randomization \cite{Gramfort12RandomizedSparsityClustering}.
%Secondly, a multi-stage procedure is proposed to gradually eliminate irrelevant voxels and reduce the sparsity penalty on the detected relevant voxels. %and multiscale point of view, which is related with random lasso that is also multi-stage \cite{Wang2011RandomLASSO}}
%In addition, we give more analysis on robustness of this algorithm by trying different clustering methods and clustering sizes, which are the key components of our algorithm.
% ; say more in details about the probability methods. }
 We make a further analysis of this ``randomized structural sparsity" based algorithm, by considering the sensitivity of several key parameters of our algorithm such as the number of the clustering, the block size of block samplings, for the first time. In addition,  we proposed a simple but effective heuristic way for setting the threshold value, which is used to filter out the uninformative features. Finally, we consider both structural MRI and functional MRI modalities.  The revealed potential biomarkers from these two modals will be put together for  better understanding of the ASD than the single modality.}

%Secondly, while we have empirically showed that %the random subsampling can also help remove the ``laziness" of multivariate classifiers in terms of discriminative pattern identification. % in  %we introduce a coarse-to-fine multi-stage procedure to gradually remove the uninterested regions to reduce the dimension of the interested voxel space  as the algorithm  moves on, in order to further accelerate the computational efficiency and enhance the stability and accuracy of support identification. %This procedure is very helpful for large scale data.
%Fifth, \comm{Where our method is better than [31] and [32], and randomized L1? Theoretically better balance of false positive and false negative}
%
%\comm{Sixth, Demonstrate
%this ``blocked" variant of stability selection can achieves significantly better sensitivity with the control of false positives for voxel selection in \cite{wang14RSS}, we show the better reproducibility and stability of selected features via our method than other alternatives based on the multi-center data. % and false negative via numerical experiments.%}
%
%TSpecifically, we perform the comprehensive .  %While the differences between the univariate and multivariate methods in terms of support identification are clearly demonstrated, the ``blocked" variant of stability selection . %Specifically,  in many existing studies, they think that due to the  strong correlations across regressors, i.e. columns of $\mathbf{X}$, univariate approaches are often more effective
%than multivariate approaches. In fact, this is not the truth, if we incorporate the idea of randomized block subsampling into the algorithm design as we do in this %paper.
%our carefully design synthetic data.  %\comm{Our methods are  exploratory, rather than confirmative.} % a wrong viewpoint.

%\item Compared with ttest
%\item Compared with pure Logistic regression
%\item Compared with SVM
%\end{itemize}


%Third,

 %lusters can be merged together ?),  rather than the regularization techniques \cite{jenatton12multiscale}.
%\begin{itemize}

%\item Why exact clustering is not necessary (for classification)?


%\end{itemize}

%For now, we focus on structured sparsity models for  either network \cite{Bellec06ClusteringNetwork} and voxel wised.
%Currently most existing methods are based on randomized LASSO and clustering. The purpose is to improve both the support identification and the prediction %scores. There are two key components, one is the clustering algorithms for grouping estimation (We often do not have prior information how the groups are %defined.), and the other one is the structured Sparsity-Inducing Penalty term design. % and the key model is similar.
%\begin{itemize}
%\item An important application is fMRI. \comm{Consider both functional homogeneity and functional connectivity}
%\item %\comm{This work could ultimately lead to a functional-connectivity-based atlas which could then be useful not only in connectivity studies but also in framing %task-based activation studies as well. There are atlases based on anatomy (MNI (Evans et al., 1992 and Talairach and Tournoux, 1988)) and cytoarchitecture %(Brodmann, 1909) but to date there is no atlas based on functional connectivity.} \cite{Shen2013Parcellation}
%\item We will show the limitation of \comm{hierarchical agglomerative clustering} in terms of losing differentiating the discriminative part and non-discriminative %part.
%\item locality is more essential than sparsity, as locality must lead to sparsity but
%not necessary vice versa
%\item Specifically, the penalization we use is
%constructed from a tree that is obtained by spatially constrained agglomerative clustering.
%\item How we do the clustering. Try to use a different clustering method.
%\item Make use of existing structured sparsity norm, or develop new ones.
%\item Supervised further reduction of the discriminative regions.
%\end{itemize}

%\begin{itemize}
%\item Sparse deep-learning algorithm for recognition and categorisation , NOT specific for functional MRI.
%\end{itemize}
%\cite{Anderson10critiqueMVPA} is critizing. We can solve its doubts. solve the  "inherent run to run variability"

%The use of subsampling  is certainly not new; we merely try to put it into a more formal framework and to show certain empirical and theoretical advantages of %doing so}
The rest of the paper is organized as follows. In section \ref{Sec:Method}, we revisit our recent algorithm for stable voxel selection. In section \ref{Sec:Experiments}, we demonstrate the advantages of our algorithm  on the real multi-center neuroimage data in terms of better consistency of selected discriminative voxels  across different centers. %both higher sensitivity and specificity.
In section \ref{Sec:conclusion}, a short summary of our work and some possible future research directions will be given. %research

\section{Randomized Structural Sparsity Applied to Multi-Center Data Analysis} \label{Sec:Method}

\subsection{The Background and Motivation }

{We focus more on how to obtain a reliable support identification in context of many fields including mental decease or external stimuli including multimedia recognition. This is not a trivial task. In this paper, we revisit  the Randomized Structural Sparsity (RSS, for shot) based algorithm we proposed in \cite{Wang14RSS}, and aim to demonstrate that this algorithm can achieve a stable and consistent feature selection result even in the existence  of data variations.

MRI data is profoundly affected by experimental and methodological factors.
There is a high likelihood of introducing undesirable inter-site variability into the data, reducing the benefit of the multi-center design. Here we take this challenge of data variation among  multiple centers as a favor to demonstrate the stability and reliability of RSS. %In particular, we consider each of the centers separately and apply the feature selection method to obtain the potential discriminative features. Then we compare the results from different sites to see whether they are consistent to certain degree. While this seems to be a natural and better way to validate the reproducibility of the final selected biomarker than only considering a  single site data,
We empirically show that  most existing feature selection methods fail to obtain a consistent result because they are not able to deal with  the possibly relatively high data variability across different centers.}

%As major advances have been made in the availability of shared neuroimaging data, such that there are more than $8,000$ shared MRI (magnetic resonance imaging) data sets available online \cite{Poldrack14neuroimaging}.


%The quantitative analysis of pooled data from related  magnetic resonance imaging  experiments has the potential to significantly accelerate progress in brain mapping.  }

%{While most of existing efforts to deal with multi-center data is to} merge the data from different sites into a pool  and analyze it as a whole, like a big single center, we are  considering the multi-center separately \cite{Martino13MP}.  Now we first review the commonly used feature selection methods in the following parts
\subsection{Review of Randomized Structural Sparsity}
%Note that we are not going to  consider how to design and implement proper protocols at different sites to reduce the variability , but
We first review the Randomized Structural Sparsity based feature selection method. %focus more on how to obtain a reliable support identification in context of many fields including mental decease or external stimuli including multimedia recognition.
%The corresponding mathematical models share the similar forms as follows. %such difficult situations. %but focus on data analysis.}
%
%\cite{Liu09:SLEP:manual}
%\comm{Interpretability: Much More than Sparsity}
%\begin{eqnarray}\label{eq:L1}
%\min_{\mathbf{w}}\|\mathbf{w}\|_1 + \lambda \|y - \mathbf{X} \mathbf{w} \|_2^2
%\end{eqnarray}
%We first review our algorithm recently proposed in \cite{Wang14RSS}. % which is expected to find out the true discriminative features  for the available high dimensional  training data with classification information in a stable way.  % in a computationally efficient and statically reliable way.
%
 %   We take the fMRI data as the example to explain our main ideas while this algorithm can also be applied to many other kinds of data such as structural MRI data. % The activation is assumed to occur in groups of voxels rather
%than in a single voxel, because the neighboring voxels within the same function area usually behavior in a similar way. In addition, the most discriminative %voxels are only a small portion of the total, and thus we need to consider the simultaneous   sparsity and clustering effects when designing algorithms. %simultaneously.  %\comm{is usually of strong correlation between neighboring voxels}.
%%
%
%
 Denote an MRI data matrix as $\mathbf{X} \in \mathbb{R}^{n\times p}$ where $n$ is the number of samples and $p$ is the number of voxels with $n \ll p$, and corresponding classification labels $\mathbf{y} \in \mathbb{R}^{n\times 1}$, where the binary classification and $\mathbf{y}_i \in \{1,-1\}$ is considered here.  % we try to find out the comprehensive set of the true discriminative voxels, i.e. this is a support identification problem.
%While our main ideas can be applied to  most existing sparsity based  multi-variate model in terms of either prediction or classification tasks,
We take the following sparse logistic regression for classification as the example  to describe the main idea of our algorithm. %show the existing difficulties and corresponding efforts. %In this section, we give a detailed description about how our algorithm is implemented. % We will try to make it clear for possible reproducible research by others.
\begin{eqnarray}\label{eq:L1}
\min_{\mathbf{w}}\|\mathbf{w}\|_1 + \lambda \sum_{i=1}^{n} \log(1+\exp(-\mathbf{y}_i(\mathbf{X}_i^{T} \mathbf{w}+c)))%\|y - \mathbf{X} \mathbf{w} \|_2^2
\end{eqnarray}
where $\mathbf{X}_i$ denotes the $i$-th row of $\mathbf{X}\in\mathbb{R}^{n\times p}$; $\mathbf{y} \in \mathbb{R}^{n\times 1}$ is the labeling information containing the classification information of each row of $\mathbf{X}$; $\mathbf{w} \in \mathbb{R}^{p\times 1}$ is the weight vector for the voxels; $c$ is the intercept (scalar). The voxels corresponding to $\mathbf{w}_i$ of large absolute value will be considered as the discriminative voxels. % \comm{We need to mention LASSO here. }
%When $p \gg n$, the plain $\ell_1$ model (\ref{eq:L1})  can identify at most only $n$ variables before it saturates \cite{Zou05ElsticNet}. % because of failing to  exploit any spatial or temporal prior information of the fMRI data. %It fails to identify clustered discriminative voxels  as expected.

%\comm{Say more about the features of our problems. and Also say the (2) is the baseline of our methods. } its rather than voxel units for activation.




%This  may not be a desirable feature for many practical problems, particularly in brain decoding.


%, i.e. we are dealing with a supervised learning problem.
%
 % Due to the %large scale of this kind problem,  efficiency is of particular importance. %algorithm.
%As mentioned above, due to the extremely high dimensionality of voxel space and usually relatively very small number of training data (samples), it is often very hard %to make reasonable inferences about the discriminative brain regions only based on the plain  $\ell_1$ based model (\ref{eq:L1}).
%though sparsity has traditionally been connected
%with interpretability, in the sense that sparser models
%are easier to interpret.
%because it  is likely to return very small number of discriminative voxels, which are not either reliable or easily readable and interpretable. Specifically, we are %focusing on voxel selection in terms of certain classification task and the common  $\ell_1$-norm based logistic regression model is as follows:
%\begin{eqnarray}\label{eq:L1}
%\min_{\mathbf{w}}\|\mathbf{w}\|_1 + \lambda \sum_{i=1}^{m} \log(1+\exp(-\mathbf{y}_i(\mathbf{w}^{T}\mathbf{X}_i+c)))%\|y - \mathbf{X} \mathbf{w} \|_2^2
%\end{eqnarray}
%where $\mathbf{X}_i^{T}$ denotes the $i$-th row of $X\in\mathbb{R}^{m\times n}$; $\mathbf{y} \in \mathbb{R}^{m\times 1}$ is the labeling information containing the %classification information of each row of $X$; $\mathbf{w} \in \mathbb{R}^{n\times 1}$ is the weight vector for the voxels; $c$ is the intercept (scalar). The voxels %correponding to $\mathbf{w}_i$ of large absolute value will be considered as the discriminative voxels. In addition, things can be even worse, due to the the strong %local correlation between columns of $\mathbf{X}$,  the model (\ref{eq:L1}) itself could  %be find the involved brain regions effectively in an interpretable way, %because
% likely makes  the plain $\ell_1$ model  select arbitrarily one of them and not the others��  resulting into  a high rate of false negatives in the
%support identification.
%There are two common ways to deal with the above situation. One is clustering and the other is randomlization based resampling or subsamping. We will give a %brief review of them in the following parts. %like to combine both in a more efficient way }

%which might not be grouped into a brain region of significant size, but scattered separately and therefore


%While the strong local correlation of columns of $\mathbf{X}$ brings difficulty for reliable estimation based on the plain $\ell_1$ model such as (\ref{eq:L1}), it also provides us with opportunity that we can identify  much more involved discriminative voxels from the only small number of samples, if we make proper use of  this correlation property, which is reflecting the  segregation and integration of the brain, i.e. the voxels are not working individually, but as groups.

% The sparsity can be  considered in terms of the groups (brain regions) of voxels, i.e. only a small fraction of all the groups are the discriminative ones, rather than in terms of the individual voxels. % and this makes it possible to
%Fortunately,
% In other words, we make our hypothesis on the spatial layout of the neural coding: sparsity, i.e. few relevant regions implied in the cognitive task, and %compact structure, i.e. relevant features grouped into connected clusters. Correspondingly we mainly focusing on the block structured sparsity, i.e. the brain %activation voxels are often clustered into several regions.
% However,  statistical analysis of fMRI data should identify cluster units rather than voxel units for activation. Cluster-based analysis was commonly applied on fMRI research to avoid the drawback that activation maps obtained by single-voxel analysis had the low SNR in individual voxel data \cite{Heller2006Cluster}. In other words,  the features in this paper are clustered voxels, which are not  exacted via  %mainly consider the problem of identification of the involved brain regions, i.e. voxel selection, and
%do not consider features extraction and their corresponding methods such as
%Principle Component Analysis \cite{Jolliffe05PCA}, Linear Discriminant Analysis \cite{Mika99FDA},  Non-negative Matrix Factorization \cite{Lee99NMF}, or  Independent Component Analysis \cite{Hyvarinen:2000:ICA,Daubechies09SparsityOvreICA}. %e  our


Structured sparsity models beyond the plain $\ell_1$ models  have
been proposed to  incorporate  more structured
constraints on the solution \cite{Bach12StructuredSparsity,Li13ClusterGuided,Mairal13SFS}. %As an important special case,   %these new structured sparsity
%models take in account the existing relationships and structures between
%the voxels (e.g., spatial, hierarchical or related to the physics of the problem at hand). %Therefore, they
%promise an even greater ease of interpretation
%of the activation patterns and much more stable results than the plain $\ell_1$ model. %, because the active voxels are
%grouped together in possibly few clusters. %, which also fits
%well with our knowledge about the brain��s specialized
%regions and networks mentioned above. That is to say, existing relationships and structures between
%the voxels (e.g., spatial, hierarchical or related to the physics of the problem at hand)
%are not disregarded, i.e. being able to encode such specific spatial constraints, potentially  for interpretability.
%and
The common way to make use of the clustering or grouping structure is the following group sparsity induced norm induced  model \cite{Bach12OSP}.
\begin{eqnarray}\label{eq:structuralL1}
\min_{\mathbf{w}}\sum_{g\in \mathcal{G}}\|\mathbf{w}_g\|_2 +  \lambda \sum_{i=1}^{n} \log(1+\exp(-\mathbf{y}_i(\mathbf{w}^{T}\mathbf{X}_i+c))),%\|y - \mathbf{X} \mathbf{w} \|_2^2
\end{eqnarray}
where $\mathcal{G}$ is the grouping information. %Compared with (\ref{eq:L1}), the main difference is the regularization term and now we are using a mixed $\ell_1/\ell_2$ norm. %\comm{Say more about the advantages}
%However,  %the group sparsity-induced norm regularized model (\ref{eq:structuralL1}) expects to improve the sensitivity  than the plain $\ell_1$ model (\ref{eq:L1}) due to the adopted mixed $\ell_{2,1}$ norm if the grouping information $\mathcal{G}$ is reliable enough.
However, the obtaining of appropriate $\mathcal{G}$ might be difficult in practice, and the final results might be too biased by the grouping information $\mathcal{G}$.   % partially because %for the following two reasons. % \comm{difficulties} the identified discriminative regions might not be accurate in terms of the given classification task. % if  the number of available training sample is small,  partially due to the very high dimensionality of the voxel space.  More important, there still exist several other problems with using of the structured induced norm.   %\comm{We do no have the structure information beforehand.}
%The main reason is that
%One reason is that
%most of methods of obtaining  $\mathcal{G}$  are not incorporating the available classification or labelling information, and  % the usually  have nothing to do with the classification task.
%it might only be derived  from either the prior anatomical knowledge or data driven methods based on the voxel correlation.
 %Moreover, due to the very high dimension and few training samples, the  finite sample control of false positives is much important, but not achieved yet.
%\comm{The other reason is that }

%, where the correlations between each activated voxel and its neighboring  activated voxels should be high while the correlations between an activated voxel and its neighboring non-activated voxels should be low. The significant discriminative regions in terms of a given classification task usually do not necessarily match the prior anatomical knowledge of brain partition, or the clustering results based on the fMRI data.

%Firstly, we usually do not have the information about the groups beforehand. Secondly, the structure or groups are only determined by unsupervised learning such as clustering based on  certain "distance" definition,  and it may not reflect the true boundaries of the involved brain regions in the given classification task.  That is to say,  t
%The grouping information may not correspond to any biological entities or certain disease, or outside stimuli, causing difficulty in interpretation.
%However, the simple clustering information based on correlation or distance might lead to probably misleading support identification results.



%The resulting penalty has
%been shown to improve the prediction performance and/or interpretability than the plain $\ell_1$ model, only
%provided that the block structure is relevant and reliable \cite{Jenatton11StructuredSparsity}, which is not necessary satisfied in practice.
% clustering is not deterministic and unique somehow, and is also usually a greedy procedure and  has high-variance.
%Due to  the underlying grouping structure is usually  %estimated by, for example,  clustering algorithms, which are greedy algorithms and be biased. Therefore, a single shot of running model (\ref{eq:structuralL1}) %might not produce satisfying results.




%While this procedure might be time-consuming, or even prohibitive for real large scale brain decoding problems, due to at least the following $3$ reason. \comm{the advantage of our algorithm over the traditional structured sparsity  is better efficiency}
%tructured sparsity where  It does not makes use of the label information. \comm{The limitation in terms of computational efficiency}
%We need averages.


%\comm{Introduce the traditionally way}
Recently, stability selection \cite{Meinshausen10StabilitySelection} based on the plain $\ell_1$ model has been widely applied \cite{Ye12SparseLearningSS,Cao2014Sparse,Ryali20123852}, due to its finite sample
%As mentioned before, an effective way to
control of the false positives.  In addition, it makes the choice of the regularization parameter insensitive.
%when applying the sparsity regularization based models is s which have been applied onto voxel selection or connection  selection .
However, it fails to make use of the prior structural information of the discriminative voxels,  and may result into a big false negative rate in order to keep low false positive rate.

Therefore, we proposed ``randomized structural sparsity" in \cite{Wang14RSS}, which  incorporates the spatial structural knowledge of voxels into the stability selection framework. It aims to  achieve  low false positive rate and false negative rate simultaneously.  In the paper, we will further reveal its stability and reliability via the multi-center data analysis, i.e. a consistent biomarker can be found across different sites due to its robustness to relatively small data perturbation. It is of important significance for the reliable scientific truth.
%
%While it has been well known that randomization allows us to handle the
%huge feature space and prevents over-fitting \cite{BangpengYao11CRD}, the
%As we have mentioned before, the other way to is to %while the adoption of the structured sparsity might produce a better result  than the traditional plain $\ell_1$ based sparse models such as LASSO,  the results of voxel selection  could still be unstable, i.e.  group structure-induced norm itself is still not enough for stable or reliable in terms of accurately finding out the comprehensive set of the involved regions, because the voxel space is of too high dimensionality and the grouping information does not take the given label information, i.e. $y$ into consideration. % \comm{For our purpose, we do not need group sparsity induced norms}
% especially for the our high dimensional cases. %ollowing several reasons. Firstly,  the underlying structure needs to be estimated by, for example, some clustering algorithms, and might be biased.
%In such cases, we will
 %make use of the ideas of
 %stability selection  is a concrete implementation in cases of variable selection, and
 % as well as actively exploring the boundaries of the involved brain regions. %For example, \cite{Varoquaux12Clustering} proposed to use randomized rescaled clustering, which means running many times of clustering %and taking a very long time. We need a computationally feasible way for large scale problems.
 %is expected to  enhance and
%improve existing variable selection methods based on sparsity such as the model (\ref{eq:L1}), and its results   %\comm{Specifically, Variable selection and structure estimation improve markedly for a %Rrange of selection
%ethods if stability selection is applied.}
%Its another advantages is that it
%are not sensitive to the sparse regularization parameter $\lambda$.  The standard stability selection  \cite{Meinshausen10StabilitySelection}  count how often each voxel is selected across the repetitions via resamping. The more often, the more likely the variables are believed to belong to the true support. \comm{A particularly attractive feature of Stability Selection is
%the error control provided by an upper bound on the
%expected number of falsely selected variables}
%
%
%Due to the fact that different variables in a strongly correlated group
%can be selected across the resampling estimates,
%Stability Selection, this is because these variables can ��split the vote��, with the result that
%neither/none of the variables is chosen}
%\comm{So we combine randomness with (\ref{eq:structuralL1})}
%
%In addition, it can remedy the the possibly incompleteness of multivariate classifiers based support identification.
%
%
%We will first briefly review the classical stability selection and point out its limitation as follows.% Tturn to making use of the grouping information to enhance the performance of stability selection in both efficiency and effectiveness. %In the case of
  %\comm{To mitigate this problem, it is possible to use randomization techniques, }
%In \cite{Gramfort12RandomizedSparsityClustering}, the authors proposed to adopt the so called randomized clustering. Specifically,
%in randomly perturbing the design matrix $\mathbf{X}$ by taking
%only a fraction of the training samples and randomly
%rescaling each variable, in our case each column of $\mathbf{X}$. Then  the clustering algorithm is run and this way
%clustering result is different for every
%randomization, as was done in the Random
%Forests algorithm \cite{Breiman01RandomForests}. And the resulted models such as (\ref{eq:structuralL1}) or its variant are also different.
%By repeating the later procedure and then counting how often each voxel is estimated to be involved, we can get a
%score for each voxel.
%\comm{The existing work} A naturel way for extending the common stability selection is that after each random rescaling of each voxel as the common stability selection, %
 % or %(\ref{eq:everaged}), as did in \cite{Varoquaux12Clustering}.
% \begin{eqnarray}\label{eq:everaged}
%\min_{\tilde{\mathbf{w}}} \|\tilde{\mathbf{w}}\|_1 +\lambda \sum_{i=1}^{m} \log(1+\exp(-y_i(\tilde{\mathbf{w}}^{T}\tilde{\mathbf{X}}_i+c)))%\|y - \mathbf{X} %\mathbf{w} \|_2^2
%\end{eqnarray}
 %
%\comm{
%Because of the randomness of the cluster shapes
%in the different iterations, the selection frequency score
%can outline the true discriminative regions, which are even   of a smaller size than the cluster
%size. %Notice that the times of rescaled clustering is the same as the repetition  times of stability selection and therefore could be large.
%} %In , it used the randomized rescaling clustering for avoiding bias and trying to outline structures not exactly matching the shapes of the clusters, for %example, those of a smaller size than the cluster
%size.
%
%  are small and connected areas  the brain  and should go along with the clustering of voxels.
%
%
%
%they assume the group information, which is usually not available beforehand. The second is the group-induced norm is not as good as it seems in this case. The %randomization bring more computational burden when clustering and group-induced norm is used. % , we are also using (averaging) as a model reduction way to
%\comm{As for the perspective of columns, or voxels} \comm{We do not need the rescaling as the original Stability selection did.}
%
%
%Now we discuss the  ``laziness" of the existing non-randomized methods in terms of  the comprehensive support identification.  It is well known that while the high classification or prediction accuracy might be achieved, most of existing multivariate classifiers including the model (\ref{eq:L1}) fail to find out all the true significant discriminative brain regions, partially because these two goals are not necessarily the same in many cases and the discriminative voxels are redundant, meaning that part of the total discriminative voxels might still achieve  not worse prediction or classification accuracy. We will demonstrate the randomization used in stability selection helps cure this kind of  ``laziness". %In other words, the current multivariate classifiers is ``active" for good prediction accuracy, but ``lazy" in recovering discriminative voxels, partially due to the existence of redundant feature sets that can lead to different models of equal predictive power.
%
%\subsection{The Key Component: Constrained Block Subsampling} \label{Sec:ConstrainedBSub}
%\comm{Randomization is an effective way of "exploration" and we re-estimate the sparse model many times perturbing the design matrix $\mathbf{X}$ or %sub-sampling data and counting how many times a given regressor (the voxel, here) is selected.  In order for better efficiency, we need to consider the %specific features of the give problems when implementing the randomization scheme. For fMRI, when subsampling or reweigthing each voxel (the columns of %$\mathbf{X}$), we need to consider the brain segregation.}
%{Need a re-writing. }
%
%\comm{The bootstrap and the other resampling methods for such dependent processes thus
%differ from their classical versions for independent and identically distributed (iid)
%variables by resampling ��blocks�� of observed variables as compared to resampling
%a single observation at a time. The intuitive reasoning behind this is that the dependence
%structure of the process is preserved within the blocks, and therefore, valid
%approximations can be obtained by resampling the blocks}
%As mentioned before, for our support identification  problems, it has several features.  The first one is that it has a very high dimensional voxel space, compared with the very limited number of training samples. For example, number of unknowns is
%commonly over hundreds of thousand of voxels, while the number of observations (treated as training data)
%never exceeds a few thousands, even hundreds. The second is that the neighboring voxels are strongly correlated reflecting the brain segregation and integration.
%Therefore, the plain $\ell_1$-norm based model does not work well in cases due to the fact that it only focuses  on cardinality and cannot easily specify or incorporate %this information about the patterns of voxels. % induced in the solution
%In this paper, we first further analyze the limitations of the two above existing solutions for this difficulty, i.e. group sparsity  based methods and randomlization %based methods. % to deal with
%Then we will propose our new algorithm which overcomes the limitations of these two kinds of existing methods and expects to obtain more stable and reliable results.
%
%
% As mentioned before, one way of obtaining more stable results are to  incorporating the grouping information into the support identification model. making use of the grouping information of voxels  provided by the clustering algorithm is to apply the group sparsity-induced norm based model like (\ref{eq:structuralL1}) as \cite{} has done. \comm{However, the drawbacks.}
%
%
% In addition, as for the sparse model, the number of estimated known is greatly reduced from the number of the voxels to the number of clusters.  Here we %emphasize that compared with  structured sparsity induced norm based sparse model, % the resulted reduced size sparse model based on averaging has better %stability in cases of very small number of training data, i.e. the rows of $\mathbf{X}$.  \comm{Why?}
%
%Clustering is for dimensional reduction and also for numerical stability due to the correlation of the data matrix $\mathbf{X}$.
%
%. % %of the brain
%reduce the size of $\ell_1$ model. \comm{Further analysis compared with group induced norms}
%\item Structured is not stable.
%\item Structured sparsity is computationally expensive for high dimensional problems.
%\end{itemize}
%In this paper, we propose an \comm{efficient} algorithm framework. \comm{in more details. the advantage of our algorithm is better efficiency}
%
%
%\begin{itemize}
%\item
%
%
%
%\comm{The second kind of methods}
% are based on stability selection \comm{} \comm{Give more detailed introduction, including rescaling}, i.e. randomization  in the context of voxel selection yields %better and stable results.
%
%
 %Specifically, we propose a variant of stability selection, which %. %computationally efficient exploring way to find out all the involved significant discriminative brain regions by making better use of the two hypothesis and overcoming the ``laziness", for large scale real fMRI data.
%Here,  we are considering both the computational challenges and the inference quality, for  the fMRI data, where are usually of extremely high dimension.
%
%Specifically, we combine  the ideas of stability selection and sparsity based voxel selection mentioned above,  together with our proposed ``constrained block subsamplings " to achieve our goals.  %\comm{Why stability selection framework?}
%
%Specifically, unlike the common stability selection where the variables or voxels are considered individually,
% takes the assumption into consideration that the discriminative voxels are distributed clusters, when performing  randomization in terms of voxels, i.e. columns of the data matrix $\mathbf{X}$. % and correspondingly
% Specifically,
%\begin{itemize}
%\item The above three difficulties, the prior . This makes exploratory models highly desirable, because there is little prior knowledge about which regions %    should be involved, \cite{NIPS2009_0802, Watanabe69Ugly-duckling-theorem}.
%    \item Clustering problems and fast convergence from the computational point of view.
%\item The third problem of ``Laziness" is popular.
%\end{itemize}
 %, \comm{via a modified stability selection incorporating the grouping features of the discriminative voxels, }%e main goal  to %While we also rely on the use of the randomization   to  find out the comprehensive set of the true involved brain regions as \cite{Gramfort12RandomizedSparsityClustering} did,   we would like to
%make our algorithm efficient and simple enough
% {Sparsity, Clustering, Comprehensive}  %, with the same goal  The goal is also. %, we turn to the exploratory way to estimate the true boundaries of the involved brain regions in terms of the available classification information $\mathbf{y}$.  We borrow the idea of stability selection and make a tailored version for our specific task, i.e. support identification based on fMRI data.
%\comm{The blocks of consecutive observations are drawn with replacement from a set of blocks.}
 %and it might still not achieve very false positive rate from limited number of resamplings, especially the grouping size is large, because of the possibly being confounded  by  uninformative variables due to chance. Even worse, the positive negate rate is likely to be  large, which can be demonstrated in the numerical experiments.
%Correspondingly,
%we propose a variant of the stability selection. %, which takes it into consideration that the discriminative voxels  usually constitute spatially distributed clusters or groups.
While ``randomized structural sparsity" is a general concept, it has a specific implementation  for voxel-wise MRI data anlaysis,  named ``Constrained Block Subsampling". % and we will review the technical insights of it. %

 ``Constrained Block Subsampling"  prefers to using the block subsampling based stability selection \cite{Beinrucker12ExtensionSS}, rather than the original reweighting based stability selection \cite{Meinshausen10StabilitySelection}. % because the former generally likely yields an
%improvement over the latter whenever the latter itself improves over
%LASSO \cite{Beinrucker12ExtensionSS}, and it  is easier for parallelization.
%
Specifically, for the training data matrix $\mathbf{X}\in\mathbb{R}^{n\times p}$, subsampling based stability selection consists in applying
the baseline, i.e. the pure $\ell_1$ regularization model such as (\ref{eq:L1}), to random submatrices of $\mathbf{X}$ of size $[n/L] \times [p/V]$, where $[]$ is  the round off to the nearest integer number, and returning those features having the largest selection frequency. The original stability selection \cite{Meinshausen10StabilitySelection} can be roughly considered as a special case of it, where $L=2$ and $V=1$, except that the original stability selection \cite{Meinshausen10StabilitySelection} reweighs each feature (voxel, here) by a random weight uniformaly sampled in $[\alpha,1]$ where $\alpha$ is a positive number, and subsampling can be intuitively seen as a more crude version
by simply dropping out randomly a large part of the
features \cite{Beinrucker12ExtensionSS}. %Authors in \cite{Beinrucker12ExtensionSS} has also showed that $L$ reflects the bias/variance tradeoff: as $L$ increases, the variance is reduced, but on the
%other hand we expect that the selection probabilities of relevant and noisy
%features become more similar and thus more difficult to separate since the
%subsample size is smaller. %However, our analysis explains in principle why it
%can be of advantage to split the sample more, provided the gain in variance
%outbalances the loss in bias.
%2.2
%
%For our problems, the relevant voxels are intercorrelated and there will be ``mutual masking", i.e. these voxels can not be simultaneously selected, if the pure $\ell_1$ model (\ref{eq:L1}) is applied.
 %While subsampling can  help solve this ``mutual masking" problem as explained in \cite{Beinrucker12ExtensionSS} to certain degree,
 %
 We further proposed to make use of the block subsampling \cite{Lahiri99Blockbootstrap},
%considering that the $\mathbf{w}$ is expected to have a sparse block structure. The block subsampling
which aims to replicate the correlation by subsampling  blocks of voxels instead of scattered voxels. %ng has been used mainly with data correlated in time (i.e. time series) but can also be used with data correlated in space, or among groups (so-called cluster data)
%\cite{Lahiri99Blockbootstrap}. %which also has extra advantages including reduced the false negatives.
%Intuitively, }
The intuition  of ``blocking" exists in the assumption that the voxels are partitioned into spatially contiguous homogeneous subgroups, though possibly distributed. %Here the ``homogeneity" %, which are ho
%is defined by the strong correlations of the voxels.

 % The intuitive reasoning behind the ``block" subsampling is that correlated voxels are assumed to be often neighbors, resulting into localized regions. %Thus the dependence or correlated structure is preserved within the blocks.
  % and this way, we will obtain the spatially contiguous discriminative regions, which matches our second hypothesis. %current knowledge about the significantly discriminative voxels. }  % by subsampling the blocks \cite{Lahiri01BlockSize}. %That is to say, the voxels are divided into relatively homogeneous subgroups or blocks.
%rather than various complicated clustering algorithms,
%for large scale problems not only due to its respecting the clustering feature of discriminative voxels, but also due to its high simplicity and efficiency.
%The simple subsampling of the classical stability selection will fail in case of fMRi data or at least not efficiently, as it is not able to replicate the %correlation in the data.
 % i.e. we use a random subset of
%voxels at a time.
%
%
%\comm{Block resampling has two purposes. Numerical stability and give an accurate shape of the identified supports.}
%In this paper, we propose to use { Block resampling is also important for the stable and comprehensive detection results}.
%The block resampling is often used in the cases of correlated data.
%In general, the block resampling is performed as follows.  One pattern of departure from independence of character evolution is correlation of evolution in adjacent characters. When this is thought to have occurred, we can correct for it by sampling, not individual characters, but blocks of adjacent characters. This is called a block bootstrap and was introduced by K\"{u}nsch (1989). If the correlations are believed to extend over some number of characters, you choose a block size, B, that is larger than this, and choose N/B blocks of size B.
%\comm{blocked sampling}
%\comm{What is our modifications to common block resampling ? Together with clustering}
%In other words, we randomly pick a small subset of all the voxels and these voxels consists of several blocks, not like the common subsampling which usually %scattered voxels.



 %A natural and naive variant of stability selection based on block subsampling can be as follows.  The total picked voxels of these   randomly decided blocks consist of only a small subset of all the voxels, for example, only $\beta p$ voxels are picked each time where $\beta=0.1$. %The block subsampling both consider the randomness in order to better explore the true boundary of the involved brain regions, and consider the spatial contagiousness of the discriminative brain regions.
%This way, we have a submatrix of the original data matrix $\mathbf{X}$, denoted as $\mathbf{X}'.$   %each random block subsampling chooses a random subset of voxels at a time. These chosen voxels are randomly distributed blocks of voxels.  %The ``laziness" of support identification of the existing multivariate classifiers is overcome because only a  subset of the total voxels are considered during each resampling.
%While there might exist many ways to make use of the available partition information $\mathcal{G}$,  we propose a simple way as follows.
 %The chosen voxels lying in a  selected random block are treated as a group $g$, and the set of all these groups  are defined as $\mathcal{G}$. Each block is a cube with length $l$, i.e. consists of $l^3$ spacial contiguous voxels for 3-D brain image.   Correspondingly, ones can solve the following group-sparsity based recovery model. {\small \begin{eqnarray}\label{eq:structuralL1Subsampling}
%\min_{\mathbf{w}'}\sum_{ g\in \mathcal{G}}\|\mathbf{w'}_{g}\|_2 +  \lambda \sum_{i \in \mathcal{J}} \log(1+\exp(-\mathbf{y}_i(\mathbf{w'}^{T}\mathbf{X}'_i+c)))%\|y - \mathbf{X} \mathbf{w} \|_2^2
%\end{eqnarray} }where $\mathbf{w}'$ and $\mathbf{X'}$ are corresponding parts of $\mathbf{w}$ and $\mathbf{X}$, respectively,  based on the selected voxels during the subsampling. % and  $\mathcal{G}$ is a predefined partitions of the brain, based on the either biological knowledge or data driven learning or estimation such as clustering.
%\comm{Say several of its advantages and disadvantages. }
%The rationality of ``blocking" exists in the assumption that the voxels are partitioned into spatially contiguous homogeneous subgroups. Here the ``homogeneity" %, which are ho
%is defined by the strong correlations of the voxels. %This way the correlation within blocks
 % and has nothing to do with the voxel correlations.
%\comm{Stability Selection is more general, not just for clustering or grouping} \comm{An important example
%is the Stability Selection methodology for variable selection proposed in Meinshausen and
%B��uhlmann (2010) and further developed in Shah and Samworth (2013). S


  %\comm{Why we think in our cases, reweighting is not as good as subsampling? be in more details. }

  %,  based on quickly exploring most probable grouping combinations of all the voxels


%as the natural supplement of sparsity in order to
%obtain interpretable inferences from almost any sparsity based supervised
%learning methods. we needs randomness \comm{in terms of both rows and columns and bring computational inefficiency}.   While the stability selection is using %the resamping techniques, which means repeating solving a sequence of sparse optimization problems
%\comm{the important role of randomness?}Randomlization is one key for stable results in the area of brain imaging. So one time
%tree construction that is obtained by spatially-constrained agglomerative clustering might not be enough \cite{jenatton12multiscale}, though this paper is more %focusing on the introduction of  a hierarchically structured regularization which takes into account the spatial and multi-scale structure of fMRI data.
%Most of the existing works ignored the resampling in terms of voxels.
%\comm{The limitation in terms of computation-- clustering costs a lot of time}


%The randomized lasso does not fully address the problem
%of large groups of correlated variables. Indeed, the
%selection score will undoubtedly tend to zero as the
%size of group increase, and be confounded by selection
%scores of uninformative variables due to chance.





%While this naive combination of stability selection and structural sparsity has the potential to perform better in cases of  large groups of correlated %variables than the original stability selection based on the plain $\ell_1$ sparse model,



%It can be viewed as using an additional 8 mm smoothing kernel on top of the spatial smoothing that was performed during preprocessing. A problem with such high level of smoothing is that the boundary of the detected clusters can no longer be trusted, especially with the liberal selection probability threshold. The tradeoffs between spatial continuity and localization should be discussed. Importantly, objective quantitative evaluations are needed.

%We are mainly focusing on the computational challenge.   %Theoretically we need to try as more as possible sets of groupings in order for stable and reliable results. %, but it is almost impossible in practice due to the very huge number of possibilities of combinatorial explosion.
%Theoretically, there have existed too many possible grouping.
%So
%If the  block subsampling version of stability selection is performed in a naive way, we might need to run too many  subsamplings in order to obtain reliable %results, and this is not affordable in practice for large scale problems. % However, from practically computational point of view, we can not afford too many.
%In order to further improve both the computational efficiency and statistical inference,
Moreover, we  incorporate the  parcelling information of  the brain into block subsampling, resulting in the so called ``constrained block subsampling'', where the ``constrained" means that the parcelling information will be respected.  The kind of partition information can be based on either the prior anatomical knowledge of brain partition \cite{TzourioMazoyer2002AnatomicalParcellation}, or the clustering results based on the MRI data. In particular, for each cluster $g \in \mathcal{G}$, it may consist of either only one or several distributed localized brain regions or called partitions, because  a cluster based on ``homogeneousness" could be disrupted into several different brain areas \cite{TzourioMazoyer2002AnatomicalParcellation}. %In terms of brain network, these areas have strong functional connections, though we do not consider brain network in this paper.  %proposed whole brain parcellation algorithm divides the brain into clusters that are homogeneous from the anatomical and functional point of view, each of which constitutes a node.
%
%
%
Correspondingly,   %each sampled block is necessarily be within a single block and it can cross different %partitions. In other word,
the selected voxels from the same cluster, after the common block subsamplings, will be considered as a subgroup, i.e.,  %from different block subsumplings. %can have
 %each random block subsampling chooses a random subset of voxels at a time. These chosen voxels are randomly distributed blocks of voxels.  %The ``laziness" of support identification of the existing multivariate classifiers is overcome because only a  subset of the total voxels are considered during each resampling.
%While there might exist many ways to make use of the available partition information $\mathcal{G}$,  we propose a simple way as follows.
 the chosen voxels lying in a  cluster $g \in \mathcal{G}$ are noted as  a set $g' \subseteq g$. Furthermore, in order to make sure that those small clusters can also be sampled during the block subsampling, we borrow some idea of ``proportionate stratified sampling" \cite{Sarndal03stratified,Devries86}, i.e. the same sampling fraction is used within each partition, in order  to reduce the false negatives, especially when the sizes of different partitions are of quite a range. %we expect that each selected block lies in a homogenous region and consists of strongly correlated voxels.
%
%If you want to be able to talk about subgroups, this may be the only way to effectively assure you'll be able to. If the subgroup is extremely small, you can use different sampling fractions (f) within the different strata to randomly over-sample the small group (although you'll then have to weight the within-group estimates using the sampling fraction whenever you want overall population estimates).
%However, we are using the block subsampling, instead of individual subsampling to make use of the spacial contiguousness. % When we use different sampling fractions in the strata, we call this disproportionate stratified random sampling. Second, stratified random sampling will generally have more statistical precision than simple random sampling.
%\comm{What kind of prior information? How we get it? How to make use of it in a better way?}
%
This way, we obtain the following group-sparsity based recovery model. %{\small
\begin{eqnarray}\label{eq:structuralL1Subsampling}
\min_{\mathbf{w}'}\sum_{g' \subset g\in \mathcal{G}}\|\mathbf{w'}_{g'}\|_2 +  \lambda \sum_{i \in \mathcal{J}} \log(1+\exp(-\mathbf{y}_i(\mathbf{w'}^{T}\mathbf{X}'_i+c)))%\|y - \mathbf{X} \mathbf{w} \|_2^2
\end{eqnarray}
%}
where $\mathbf{w}'$ and $\mathbf{X'}$ are corresponding parts of $\mathbf{w}$ and $\mathbf{X}$, respectively,  based on the selected voxels during the subsampling, and  $\mathcal{G}$ is a predefined or estimated partitions of the brain. $\mathcal{J}$ is the set of the indices of the selected samples of the current subsamling.

While  ``constrained block subsamplings"  respects the prior knowledge $\mathcal{G}$, it also provides the flexibility that discovered discriminative regions can be of any shape. The final selected voxels from a cluster can be only part of  it, because   the randomness of the block subsampling in the different iterations of the stability selection procedure, makes the selection frequency score
 be able to outline structures of the true discriminative regions.  This kind of flexibility is of importance because  the group information  might not exactly reflect the true shapes of the discriminative brain regions. The grouping information  is usually obtained via certain clustering algorithms which often only depend on the correlation information of voxels.
%can be part of each clusters.
% In addition, the adoption of the prior knowledge $\mathcal{G}$ helps partially preserve the boundary of the detected clusters for the tradeoffs between spatial continuity and localization. %are not necessarily matching the  the and efficiency of random exploration of all the possible groupings.  It can also be considered as a specific randomized structural sparsity, which aims to solve the possible bias of estimated grouping information.   %Notice that performing unsupervised clustering based on the strong local correlation of neighboring voxels of fMRI data is just one means to get this partition information and we take it as an example in this paper.

%\comm{Compared with results of one-shot structural sparsity model,  the idea of resamplings or randomness overcomes the heavy dependency of its results of support identification on the prior knowledge of grouping $\mathcal{G}$, which might be of big bias compared with the  mapping of the true significant discriminative voxels in terms of the given specific classification task.
%}

%The above is a general procedure and can be considered as stability selection applied to the structured sparsity based various classifiers such as SVM or %logistic regression and so on.  This constrained blocked subsampling is a balance between flexibility and efficiency. i.e  it is doing random exploration based %on the prior information about the partition.
%\comm{Say more about the one time partition is enough. Why?} \comm{Define } Subsampling allows certain estimation errors when defining the partitions $\mathcal{G}$.

%In addition, the subsampling can help overcome the
%%\comm{However,
%the Elastic Net and Structural Sparsity, which \comm{aims to reduce false negatives},  share the common difficulty in terms of dealing with the second hypothesis, i.e.  obtaining the true grouping information of the discriminative voxels. For Elastic Net, the additional $\ell_2$ term is \comm{to avoid redundancy.} only encouraging   a grouping effect, but far away from obtaining the accurate clustering information of these discriminative vovels.  %In addition,, where strongly
%correlated predictors tend to be in or out of the model together}
 %For structural sparsity,
% the grouping information is usually obtained via certain clustering algorithms and they are only unsupervised procedures without using the classification information of the training data, and only depending on the correlation information of voxels and therefore this kind of clustering or grouping result might be slightly far away from that of the true discriminative regions. %In addition, the clustering result often is of high-variance since most clustering algorithms based on greedy approaches to non-convex problems. %tree constriction, %not supervised %procedure and the label information of the training samples is not used, the clustering %or the grouping information may not be of enough %accuracy .%}
% }

%\end{itemize}
%  \begin{eqnarray}\label{eq:everaged}
%\min_{\tilde{\mathbf{w}}} \|\tilde{\mathbf{w}}\|_1 +\lambda \sum_{i=1}^{m} \log(1+\exp(-y_i(\tilde{\mathbf{w}}^{T}\tilde{\mathbf{X}}_i+c)))%\|y - \mathbf{X} %\mathbf{w} \|_2^2
%\end{eqnarray}

%Here, clustering is not only for for numerical stability, but also for dimensional reduction.
%Due to the number of clusters $p$ being much smaller than the dimension of the voxel space $n$, we have a much better chance to obtain a stable result.
%Notice that it is reasonable to adopting averaging, instead of group structure induced norm as the following model (\ref{eq:structuralL1Subsampling}), because

%we are only caring about the positions  are the components of the $\mathbf{w}$ of large magnitudes, rather than their specific values. In other word, we are %only caring about support identification in this paper. The advantages of (\ref{eq:everaged}) over (\ref{eq:structuralL1Subsampling}) is the lower dimension %and more stable results while  (\ref{eq:structuralL1Subsampling}) is expected to give more detailed results than (\ref{eq:everaged}). %Another advantage of the
%averaging method is   it further increases the stability of the support identification results�� \comm{because}. %though it indeed greatly reduces the correlation of the data mat
 %We do it together with a one-short clustering. \comm{New block resampling}We need to take the spatial correction into consideration in order for the improved stability of the sparse methods. \comm{We need the exploration, because strong spacial correlation only
%captures pairwise information and provide with rough clustering information. However, }
%\begin{eqnarray}\label{eq:blocksampling}
%\min_{\mathbf{w}}\sum_{g' \subset g\in \mathcal{G}}\|\mathbf{w}_{g'}\|_2 +  \lambda \sum_{i=1}^{m} \log(1+\exp(-y_i(\mathbf{w}^{T}\mathbf{X'}_i+c)))%\|y - \mathbf{X} %\mathbf{w} \|_2^2
%\end{eqnarray}


%\end{itemize}

%\comm{The below is all about the second aspect of block sampling}
%\begin{itemize}
%\item \comm{Plus a complimentary idea of ISD} Greedy is GOOD. Different way to use greedy ideas.
%\item \comm{
%By making use of the block resampling, the advantage of our algorithm is to implicitly and automatically take advantage of the locally strong spatial correlation and spatial multivariate patterns of fMRI. % random way, % in a %sparsity when the sparsity structure is not known in advance, to build accurate models from relatively small sample sets.
%
%As mentioned before, the success of our method is the blocked resampling:
%With
%high probability, the blocks  selected during each iteration of the randomization (stability selection) procedure
%all contain the true discriminative voxels. When backprojecting
%the selection frequency in the non-clustered
%voxel space, % than the cluster
%size.
%In addition, the number of inferred discriminative voxel clusters are not necessarily equal to the cardinality of

%its another key advantage is the computational efficiency where the number of resamplings are not necessarily very large as demonstrated in the numerical experiments.
%this block random subsumpling does not consider the prior anatomical knowledge of brain partition (\cite{TzourioMazoyer2002AnatomicalParcellation}), or the clustering results based on the fMRI data, and thus some blocks might cross different homogeneous brain regions.



%In short, ``constrained block subsampling" is the combination of structure sparsity and stability selection, where the structure information is mostly reflected in $\mathcal{G}$, which is very important for a better result \cite{Yu2012259,Huang:2009}. %In addition, this clustering information helps to reduce the possible combinations of nonzeros of $\mathbf{w}$ for subsamplings and provide an efficient  exploration of the solution space.  Moreover, considering the spacial contiguousness structure of discriminative voxels, the block subsampling is preferred instead of traditional individual sampling.   A similar idea has been used in the terms of identifying differences in brain networks based on the mass-univariate testings \cite{Zalesky2010NBS} via cluster-based statistical methods considering that connections comprising the contrast or effect of interest are interconnected. It has showed its advantages by identifying an expansive dysconnected subnetwork which a mass-univariate analysis fails to find for schizophrenia.
%Therefore, a more comprehensive set of discriminitive voxels is expected to be revealed than the alternatives with satisfying accuracy requirement.

For the specific multi-center data analysis, the inherent variability between different centers \cite{Friedman2006varability,Friedman06JMRI} requires our algorithm to pay more attention to the reduction of variance of feature selection results, though it may result in certain  bias increasing, due to the bias-variance dilemma or bias-variance tradeoff  \cite{Geman:1992:NNB}.
%
%For our case of small samples and very high dimensional feature space, we need to consider .
In general, we would like to pay a little bias to save a lot of variance. Considering dimensionality reduction can decrease variance by simplifying models \cite{James13StasiticalLearning},   %while we can still use the (\ref{eq:structuralL1Subsampling}) as the baseline subproblem for our stability selection framework,
%Randomization has been widely accepted as a key in solving high dimensional problems, not just limited to our specific problems.
%Now it comes to the ``constrained block subsamplings" for reliable support identification.
%Denote the number of resamplings as $K$.
%This blocked variant of stability selection is the same as the classical stability selection in terms of dealing with the samples, or the rows of the data  %matrix $\mathbf{X}$ by taking
%only a fraction of the training samples.
we still use the   ``averaging" idea \cite{Varoquaux12Clustering} applied to  (\ref{eq:structuralL1Subsampling}), % in regression for the fMRI data, which poses the challenge of having far more features than samples.
because \cite{Park01042007} has proved that when the variables or features were positively correlated, their average was a strong feature, yielding a fit with lower variance than the individual variables. Specifically,
%
%Let the sub-sampling fraction be $\alpha \in [0,1]$ and the indices of selected rows are denoted as $\mathcal{J}$ and the cardinality of $\mathcal{J}$ is $[\alpha n]$, where $[]$ is  the round off to the nearest integer number.
% and randomly
%picking a fraction of all the voxels.
%
%The ``block subsamplings" %refers to the resampling of the voxels, i.e. the columns of $\mathbf{X}$.
%Considering that the voxels are strongly locally correlated or the involved voxels are often grouped together due to the mutual interactions, we use the block subsampling in this randomization procedure.
%We
%randomly picks several blocks of voxels, instead of individual voxels as the common subsampling does.  The locations of
%the blocks are randomly decided and the total picked voxles consist of only a small subset of all the voxels, for example, only $\beta p$ voxels are picked %each time where $\beta=0.1$. %The block subsampling both consider the randomness in order to better explore the true boundary of the involved brain regions, and consider the spatial contagiousness of the discriminative brain regions.
%This way, we have a submatrix of the original data matrix $\mathbf{X}$, denoted as $\mathbf{X}'.$  Then the estimated parcelling information $\mathcal{G}$ is %incorporated and we have  %Notice that while  the voxels  corresponding to the columns of $\mathbf{X}$ here) are reweigthed by a random weight uniformly sampled in $[\gamma,1]$ $(\gamma \ge 0)$ in the traditional stability selection, we are in fact using the binary weights (either $0$ or $1$) here via subsamplings.
%
%That is to say, our
%\comm{In this paper, we use a binary weight, i.e. the weight is either $0$ or $1$. }
%
%``Constrained Block Subsampling" refers to  how to define the structure when applying a structural sparsity model. It makes use of both the partition information and the subsampling to define the structure as mentioned in Section \ref{Sec:ConstrainedBSub}. The model  \eqref{eq:structuralL1Subsampling} is a simple way our proposed, which
%For each block subsampling, we can solve
%is a group-sparsity based recovery model. % brings flexibility for outlining the
 %structures of true discriminative regions,
by averaging the voxels picked by the block subsampling lying in the same cluster as a single super-voxel,  the model  (\ref{eq:structuralL1Subsampling}),  can be further reduced to the following reduced dimensional version
\begin{eqnarray}\label{eq:everaged}
\min_{\tilde{\mathbf{w}}}\sum_{g' \subset g\in \mathcal{G}}|\tilde{\mathbf{w}}_{g'}| +  \lambda \sum_{i \in \mathcal{J}} \log(1+\exp(-y_i(\tilde{\mathbf{w}}^{T}\mathbf{\tilde{X}}_i+c)))%\|y - \mathbf{X} \mathbf{w} \|_2^2
\end{eqnarray}
where $\tilde{\mathbf{w}} \in \mathbb{R}^q$, and $q$ is the number of clusters. $\tilde{\mathbf{w}}_{g'} $ is an average of voxels in the subset $g'$ of cluster $g\in \mathcal{G}$, and $\mathbf{\tilde{X}}\in \mathbb{R}^{[\alpha n] \times p}$ is the corresponding averaged $\mathbf{X}$.
%}
  %The identification process can be accelerated
 % , \comm{because} we are only caring about the position information  of  the components of the $\mathbf{w}$ of large magnitudes, rather than their specific values.
% Thus the number of variables of the sparse recovery model (\ref{eq:L1}) is greatly reduced to the number of clusters and the properties of the resulting new data matrix $\tilde{\mathbf{X}}$ is also improved due to de-correlation via the clustering. %The idea of averaging was also applied in \cite{Gramfort12RandomizedSparsityClustering}. %, for purpose of obtaining the even more stable and reliable results resulting from the reduction of  the dimensional of the problem. Moreover,  the  lower dimension leads into a much faster computation. % In fact, the averaging is acting like linear combinations of original voxels, as the PCA does. % The advantages of (\ref{eq:everaged}) over (\ref{eq:structuralL1Subsampling}) is  and more stable results while  (\ref{eq:structuralL1Subsampling}) is expected to give more detailed results than (\ref{eq:everaged}). %Another advantage of the  % we prefer using the averaging of the clustered voxels to get a reduced size $\ell_1$ model as \cite{Varoquaux12Clustering} did, in this paper.
%The reduced formula has been given in (\ref{eq:everaged}). % is as follows. \begin{eqnarray}\label{eq:everaged1}
%\min_{\tilde{\mathbf{w}}}\sum_{g' \subset g\in \mathcal{G}}|\tilde{\mathbf{w}}_{g'}| +  \lambda \sum_{i \in \mathcal{J}} %\log(1+\exp(-y_i(\tilde{\mathbf{w}}^{T}\mathbf{\tilde{X}}_i+c)))%\|y - \mathbf{X} \mathbf{w} \|_2^2
%\end{eqnarray}
%where $\tilde{\mathbf{w}} \in \mathbb{R}^p$, and $p$ is the number of clusters. $\tilde{\mathbf{w}}_{g'} $ is an average of voxels in the subset $g'$ of group $g\in %\mathcal{G}$, and $\mathbf{\tilde{X}}\in \mathbb{R}^{[\alpha n] \times p}$ is the corresponding averaged $\mathbf{X}$.
%\comm{%In order to accelerate the computation, we can use the averaging idea.  We give a rough introduction here and the details will be presented in section \ref{Sec:AlgIntro}.  Averaging is like PCA, for dimensional reduction for stable results.
%\comm{Our blocked stability selection is general}
%\begin{itemize}
%\item
 %   The reduced problem  is as follows.
%
  %\comm{
%  Notice that the ``averaging" of sumsampling is more than a simple spatial smoothing, due to different sumsampling results of different stability selection iterations. Therefore, the boundaries of the detected discriminative regions can be trusted to certain accuracy. %The tradeoffs between spatial continuity and localization should be discussed. Importantly, objective quantitative evaluations are needed. %Also, comparisons against [31] and [32] as well as existing techniques that account for spatial continuity are vital.}
%\comm{Say more about our method in an even bigger background}
%
If the $j$-th column of $\tilde{\mathbf{X}}$ is selected due to the large magnitude of $\tilde{\mathbf{w}}_j$, then its represented picked blocked voxels lying in the group $g^{(j)} \in \mathcal{G}$ ($j=1,2,\ldots, q$) of $\mathbf{X}$ are all counted to be selected, in % Back projecting with
%high probability, clusters selected during each iteration
%all contain an original non zero feature. When backprojecting
%the selection frequency in
the non-clustered
space, and its corresponding score $\mathbf{s}_i$ will be updated ($i=1,2,\ldots,p$). Notice that the ``averaging" of sumsampling is more than a simple spatial smoothing, due to different sumsampling results of different stability selection iterations. Therefore, the boundaries of the detected discriminative regions can be still trusted to certain accuracy. %The tradeoffs between spatial continuity and

%In summary, we propose to use the variant of stability selection based on constrained blocked sampling together with multi-stage refineness scheme for fast support identification, especially for large scale (high dimensional) real problems.


%\comm{To avoid the possible bias brought by the one shoot of clustering, %d obtain a more stable and reliable identification result based on the pure $\ell_1$ %model (\ref{}),
%ones can turn to the help of randomization}

%\comm{Common staility is random reweighting. Ours is subsamplings. }
%\comm{The term ��bootstrap �� was coined by Efron (1979). He described both the nonparametric and parametric bootstrap. In particular, his nonparametric bootstrap is the procedure of resampling with replacement from the original sample at the same sample size, which is by far the most commonly used bootstrap procedure. It wasn��t long before people experimented with resampling at different sample sizes. But the key discovery in that area came later. Politis and Romano (1994) described resampling without replacement from the original sample at smaller than the original sample size. This is different enough from Efron��s idea that in their book (Politis, Romano, and Wolf, 1999) they don��t call it ��bootstrap �� but just plain ��subsampling��. Whatever you call it, here��s why it is such an important innovation. ? Politis and Romano��s subsampling bootstrap takes samples without replacement of size b from the original sample of size n, generally with b <<n (read ��b much less than n��). Such samples are themselves}


%\comm{Averaging and Subsampling bring advantages for the high dimensional and small samples cases.}
%\subsection{Algorithmic framework} \label{Sec:AlgIntro}

%{In summary, The parcellation is kind of bias using reduce the variance. bias variance tradeoff }

% This section  gives an overall introduction to the algorithmic framework, which can be a multi-stage scheme and the number of stages is denoted as $N$. For each stage, it consists of three parts, i.e. ``obtaining partitioning information", ``running stability selection based on constrained block subsampling'', and ``removing the least important partitions", respectively. %  which will be described in  details in the following parts.

%While the  partitioning information of voxels might be obtained through the prior anatomical knowledge, we will take the data-driven method in this paper, though their combination might be a better choice.
%Specifically, we first perform a clustering operation aiming to partition the voxels into many patches according to their strong local correlations.% While we do not prefer any  specific algorithm and  ones can adopt k-means \cite{Hartigan79Kmeans} and hierarchical clustering algorithms \cite{Murtagh85},  %, but the general task to be solved. It can be achieved by various algorithms that differ significantly in their notion of what constitutes a cluster and how to efficiently find them.  Here we do not specify any kind of clustering algorithm,  and ones can adopt k-means \cite{Hartigan79Kmeans} and hierarchical clustering algorithms \cite{}, though particular
%it might be better that the clustering algorithm will take spacial contiguousness as well as other available brain segregation knowledge into consideration.
%We denote the set of the groups via the clustering algorithm as $\mathcal{G}$, whose cardinality is denoted as $q$, which is usually much less than $n$ and comparable with $m$. %\comm{ Since the clustering could be erroneous and highly depends on the number of clusters, introducing variability through perturbing the clustering as in [31] and [32] could be advantageous. We need to say the difference between our method with theirs} %The clustering provides us with the voxel mergering or voxel averaging within the same cluster and reduce the numerical difficulty brought by the strong local correlation of $\mathbf{X}$ in the later computation.
%In addition, the clustering information helps filter out more important groupings for the variant of stability selection based on block subsampling, in order %to quickly obtain reliable and stable results, especially when $N\times K$ is not very large. %Let us revisit the above algorithm and we would like to give more analysis on several key components and try to extend our idea into a more general setting.
%\begin{itemize}
%\item

%about the choice of clustering algorithms, %\comm{Introduce the procedure of our algorithm in more details.}
%As we know, cluster analysis itself is not one specific algorithm, but the general task to be solved. It can be achieved by various algorithms that differ significantly %in their notion of what constitutes a cluster and how to efficiently find them.  Here we do not specify any particular kind of clustering algorithm,   though it might %be better that the clustering algorithm will take spacial contiguousness as well as other available brain segregation knowledge into consideration.


 %Randomization has been widely accepted as a key in solving high dimensional problems, not just limited to our specific problems.
%Then it comes to the ``constrained block subsamplings" for reliable support identification, and  denote the number of resamplings as $K$.
%This blocked variant of stability selection is the same as the classical stability selection in terms of dealing with the samples, or the rows of the data  matrix $\mathbf{X}$ by taking
%only a fraction of the total training samples.
%Let the sub-sampling fraction be $\alpha \in [0,1]$ and the indices of selected rows are denoted as $\mathcal{J}$ and the cardinality of $\mathcal{J}$ is $[\alpha n]$, where $[]$ is  the round off to the nearest integer number.
% and randomly
%picking a fraction of all the voxels.
%
%Then ``contrained block subsamplings" are applied to the voxels, i.e. the columns of  $\mathbf{X}$ as mention in last section. %refers to the resampling of the voxels, i.e. the columns of $\mathbf{X}$.
%Considering that the voxels are strongly locally correlated or the involved voxels are often grouped together due to the mutual interactions, we use the block subsampling in this randomization procedure.
%We
%randomly picks several blocks of voxels, instead of individual voxels as the common subsampling does.  The locations of
%the blocks are randomly decided and the total picked voxles consist of only a small subset of all the voxels, for example, only $\beta p$ voxels are picked each time where $\beta=0.1$. %The block subsampling both consider the randomness in order to better explore the true boundary of the involved brain regions, and %consider the spatial contagiousness of the discriminative brain regions.
%This way, we have a submatrix of the original data matrix $\mathbf{X}$, denoted as $\mathbf{X}'.$
%hen the estimated parcelling information $\mathcal{G}$ is incorporated and we have  %Notice that while  the voxels  corresponding to the columns of $\mathbf{X}$ here) are reweigthed by a random weight uniformly sampled in $[\gamma,1]$ $(\gamma \ge 0)$ in the traditional stability selection, we are in fact using the binary weights (either $0$ or $1$) here via subsamplings.
%
%That is to say, our
%\comm{In this paper, we use a binary weight, i.e. the weight is either $0$ or $1$. }
%
%``Constrained Block Subsampling" refers to  how to define the structure when applying a structural sparsity model. It makes use of both the partition information and the subsampling to define the structure as mentioned in Section \ref{Sec:ConstrainedBSub}. The model  \eqref{eq:structuralL1Subsampling} is a simple way our proposed, which
%For each block subsampling, we can solve
%is a group-sparsity based recovery model. % brings flexibility for outlining the
 %structures of true discriminative regions,
 %the model  (\ref{eq:structuralL1Subsampling}), which  can be further reduced
  %The identification process can be accelerated
 % by averaged the voxels picked by the block subsampling lying in the same group as a single voxel, \comm{because} we are only caring about the position information  of  the components of the $\mathbf{w}$ of large magnitudes, rather than their specific values. Thus the number of variables of the sparse recovery model (\ref{eq:L1}) is greatly reduced to the number of clusters and the properties of the resulting new data matrix $\tilde{\mathbf{X}}$ is also improved due to de-correlation via the clustering. The idea of averaging was also applied in \cite{Gramfort12RandomizedSparsityClustering}. %, for purpose of obtaining the even more stable and reliable results resulting from the reduction of  the dimensional of the problem. Moreover,  the  lower dimension leads into a much faster computation. % In fact, the averaging is acting like linear combinations of original voxels, as the PCA does. % The advantages of (\ref{eq:everaged}) over (\ref{eq:structuralL1Subsampling}) is  and more stable results while  (\ref{eq:structuralL1Subsampling}) is expected to give more detailed results than (\ref{eq:everaged}). %Another advantage of the  % we prefer using the averaging of the clustered voxels to get a reduced size $\ell_1$ model as \cite{Varoquaux12Clustering} did, in this paper.
%The reduced formula has been given in (\ref{eq:everaged}). % is as follows. \begin{eqnarray}\label{eq:everaged1}
%\min_{\tilde{\mathbf{w}}}\sum_{g' \subset g\in \mathcal{G}}|\tilde{\mathbf{w}}_{g'}| +  \lambda \sum_{i \in \mathcal{J}} %\log(1+\exp(-y_i(\tilde{\mathbf{w}}^{T}\mathbf{\tilde{X}}_i+c)))%\|y - \mathbf{X} \mathbf{w} \|_2^2
%\end{eqnarray}
%where $\tilde{\mathbf{w}} \in \mathbb{R}^p$, and $p$ is the number of clusters. $\tilde{\mathbf{w}}_{g'} $ is an average of voxels in the subset $g'$ of group $g\in %\mathcal{G}$, and $\mathbf{\tilde{X}}\in \mathbb{R}^{[\alpha n] \times p}$ is the corresponding averaged $\mathbf{X}$.
%\comm{%In order to accelerate the computation, we can use the averaging idea.  We give a rough introduction here and the details will be presented in section \ref{Sec:AlgIntro}.  Averaging is like PCA, for dimensional reduction for stable results.
%\comm{Our blocked stability selection is general}
%\begin{itemize}
%\item
%    The reduced problem  is as follows. \begin{eqnarray}\label{eq:everaged}
%\min_{\tilde{\mathbf{w}}}\sum_{g' \subset g\in \mathcal{G}}|\tilde{\mathbf{w}}_{g'}| +  \lambda \sum_{i \in \mathcal{J}} %\log(1+\exp(-y_i(\tilde{\mathbf{w}}^{T}\mathbf{\tilde{X}}_i+c)))%\|y - \mathbf{X} \mathbf{w} \|_2^2
%\end{eqnarray}
%where $\tilde{\mathbf{w}} \in \mathbb{R}^p$, and $p$ is the number of clusters. $\tilde{\mathbf{w}}_{g'} $ is an average of voxels in the subset $g'$ of group %$g\in \mathcal{G}$, and $\mathbf{\tilde{X}}\in \mathbb{R}^{[\alpha n] \times p}$ is the corresponding averaged $\mathbf{X}$.
%}

%The above subsampling is followed by an averaging procedure. The picked voxels lying in the same group $ g^{(j)}  \in \mathcal{G}$ ($j=1,2 \ldots, q$) is averaged, i.e.
%the corresponding columns of the data matrix $\mathbf{X}$ is averaged and merged into one single column; the unpicked columns of $\mathbf{X}$ is
%removed from the new data matrix, and this way we have a new data matrix $\mathbf{X}'$, which has a much better conditioning and de-correlation property.

%In addition, the resulted $\ell_1$ model (\ref{eq:everaged}) has a much smaller size due to $q \ll p $ and can be solved very quickly. This point is very important for stability selection, because we need to repeatedly solve such kind of problems many times.
%Estimate $\tilde{\mathbf{w}} \in \mathbb{R}^q$ from $\tilde{X}$ and $\tilde{y}$ with the sparse logistic regression model (\ref{eq:everaged}).


%The above is the typical steps involved in the modified stability selection procedure. By
%repeating thses steps and then counting how
%often each voxel is selected across the repetitions,
%each voxel can be assigned a score. Higher scores
%denote voxels likely to belong to the true support.

%Then
%While our proposed method has already achieves a good computational efficiency, ones can adopt the following greedy multi-stage procedure as follows, in order
%``removing the least important partitions" is a greedy option and can be performed  as follows.  %  A similar idea has already been proposed by us in %\cite{Wang10ISD}.
%After one stage of stability selection, the algorithm check the average value of $\mathbf{s}_i$ for each $g^{(j)} \in \mathcal{G}$ ($j=1,2,\ldots, q$). If certain voxels are considered to be probably  not belonging to the true support and  we will discard it in the next stage in order to have a reduced data matrix $\mathbf{X}$. In the next stage, we can run the clustering again on the remaining voxels and repeat the randomization, and then remove more irrelevant voxels, for better refineness. %We can still prefer the same number of clusters as the previous stage.
%Due to the discarding of parts of voxels, we expect to obtain smaller sizes of clusters, and we can correspondingly reduce the block size when subsampling.
%These procedure is repeated and consists of a multi-stage coarse-to-fine greedy procedure, especially suitable for extremely  high dimensional problems. % The full procedure of our algorithm is plotted in Figure \ref{Fig:procedure}. %  } %\comm{It shares  some similarity with Random
%LASSO \cite{Wang2011RandomLASSO} where LASSO is used as a baseline, especially in its first phase.  In the first phase, LASSO is run repeatedly on bootstrap samples
%of the data and random subsets of the features. But for random LASSO, in a second phase, LASSO is run
%again repeatedly but this time features are drawn randomly with probability
%proportional to their importance as estimated in the first phase.}


%, this multi-state procedure can help greatly reduce the times of randomization, required for exploring and identifying the true support
%In our numerical example, we did not use multi-stage scheme, because the one-stage implementation ($N=1$) of our algorithm has already achieves an impressive %results, and we are mainly focus on the explanation of the constrained subsampling in terms of efficient support identification.
%\comm{At the beginning, all the voxels are divided into $q$ clusters.  Traditionally it is followed by the randomization based stability selection, where the (\ref{eq:everaged}) or (\ref{eq:structuralL1Subsampling}) needs to be run $\ell$ times (for example, $\ell=1000$) in order to shape out the true involved brain regions. If  the multi-stage procedure is adopted, after $K$ times randomization (for example, $K=20$), we have the rough estimate about which clusters are irrelevant.  Then we do not consider these voxels in the next stage. This way, the size of the problem is reduced and the set of possible groups of voxels is also reduced.  .
%}


%In the next stage, we perform the clustering again on the remaining voxels, on the reduced new data matrix, still denoted as $\mathbf{X}$ for simplicity. %This is optional.
 %In summary, our multi-stage procedure is a coarse-to-fine procedure taking both accuracy and efficiency into consideration.




%because of the randomness of the cluster shapes
%in the different iterations, the selection frequency score
%can outline structures of a smaller size than the cluster
%size.

% An example of a floating figure using the graphicx package.
% Note that \label must occur AFTER (or within) \caption.
% For figures, \caption should occur after the \includegraphics.
% Note that IEEEtran v1.7 and later has special internal code that
% is designed to preserve the operation of \label within \caption
% even when the captionsoff option is in effect. However, because
% of issues like this, it may be the safest practice to put all your
% \label just after \caption rather than within \caption{}.
%
% Reminder: the "draftcls" or "draftclsnofoot", not "draft", class
% option should be used if it is desired that the figures are to be
% displayed while in draft mode.
%
%\begin{figure}[!t]
%\centering
%\includegraphics[width=2.5in]{myfigure}
% where an .eps filename suffix will be assumed under latex,
% and a .pdf suffix will be assumed for pdflatex; or what has been declared
% via \DeclareGraphicsExtensions.
%\caption{Simulation Results.}
%\label{fig_sim}
%\end{figure}

% Note that IEEE typically puts floats only at the top, even when this
% results in a large percentage of a column being occupied by floats.


% An example of a double column floating figure using two subfigures.
% (The subfig.sty package must be loaded for this to work.)
% The subfigure \label commands are set within each subfloat command,
% and the \label for the overall figure must come after \caption.
% \hfil is used as a separator to get equal spacing.
% Watch out that the combined width of all the subfigures on a
% line do not exceed the text width or a line break will occur.
%
%\begin{figure*}[!t]
%\centering
%\subfloat[Case I]{\includegraphics[width=2.5in]{box}%
%\label{fig_first_case}}
%\hfil
%\subfloat[Case II]{\includegraphics[width=2.5in]{box}%
%\label{fig_second_case}}
%\caption{Simulation results.}
%\label{fig_sim}
%\end{figure*}
%
% Note that often IEEE papers with subfigures do not employ subfigure
% captions (using the optional argument to \subfloat[]), but instead will
% reference/describe all of them (a), (b), etc., within the main caption.


% An example of a floating table. Note that, for IEEE style tables, the
% \caption command should come BEFORE the table. Table text will default to
% \footnotesize as IEEE normally uses this smaller font for tables.
% The \label must come after \caption as always.
%
%\begin{table}[!t]
%% increase table row spacing, adjust to taste
%\renewcommand{\arraystretch}{1.3}
% if using array.sty, it might be a good idea to tweak the value of
% \extrarowheight as needed to properly center the text within the cells
%\caption{An Example of a Table}
%\label{table_example}
%\centering
%% Some packages, such as MDW tools, offer better commands for making tables
%% than the plain LaTeX2e tabular which is used here.
%\begin{tabular}{|c||c|}
%\hline
%One & Two\\
%\hline
%Three & Four\\
%\hline
%\end{tabular}
%\end{table}


% Note that IEEE does not put floats in the very first column - or typically
% anywhere on the first page for that matter. Also, in-text middle ("here")
% positioning is not used. Most IEEE journals use top floats exclusively.
% Note that, LaTeX2e, unlike IEEE journals, places footnotes above bottom
% floats. This can be corrected via the \fnbelowfloat command of the
% stfloats package.



%\subsection{Revisiting and Analysis of the Algorithm }
%For the fMRI data, voxels might be strongly locally correlated. Without considering this correlation, the voxel selection based on the plain $\ell_1$ model (\ref{eq:L1}) might produce kind of arbitrary and very unstable result, which is not interpretable.  In order to reduce this negative effect caused by correlation, ones proposed to perform the clustering first and get the grouping information to perform the group sparsity model (\ref{eq:structuralL1}). \comm{However, .}  Due to the very high voxel space, we propose to make use of stability selection to further improve the quality of voxel selection. The stability Selection is based on the randomized subsamping, in terms of both samples (rows of $\mathbf{X}$) and voxels (columns of $\mathbf{X}$).

%\item
%While the procedure of our algorithm is summarized in the Table of section \ref{Sec:AlgIntro}, it can also showed in gbe plotted in Figure \ref{Fig:procedure}. % in the
%\begin{figure}[h!]
%    \centering
%    \includegraphics[width=0.5\textwidth]{procedure.eps}
%    \caption{The Graphic Procedure of Our Algorithm: A variant of stability selection is proposed via blocked sampling and the incorporation of the estimated grouping information. It parters with a multi-stage refineness scheme where  for the irrelevant groups of voxels are eliminated. The detailed description of the procedure will be presented in Section \ref{Sec:AlgIntro} }
%    \label{Fig:procedure}
%\end{figure}
%next section.

%In summary,  we emphasize the role of randomization used in stability selection, which is the key for us to accurately find out all the involved brain regions.  The randomization is performed via the constrained block subsampling, which make full use of the mentioned two hypothesis above.  There have existed some arguments that complex multivariate  methods generally do not necessarily  outperform simple univariate feature selection methods \cite{Haury11biomarker}. We will show that by use of randomness, we can make multivariate methods no longer "lazy" and have much better chance  to find out the comprehensive set of the involved brain regions and behavior much better than the univariate methods, which fails for multi-voxel discriminative patterns. In addition, by considering the features of the fMRI data, our algorithm is expected to achieve  much better computational efficiency and inference quality. %compared with the common stability selection.

 %in our cases.
%This point  constitutes one of our novelties, though mostly in the empirical points of view.  % \comm{The theoretical verification is based on the previous study. about stability selection. The limitation on the ``laziness" of classifiers.}

 %Like the common stability selection is compatible with most sparse models, our variant of stability selection is  compatible with most of existing structural sparsity models. % } %  The size of the block depends on the specific problems and we find the our final results are not very sensitive to the size of blocks. %In addition, in order to accelerate the
%\item emphasize the generalization of our idea
%\end{itemize}

%Our novelty is mainly on the randomized subsampling of voxels.  \comm{We propose an efficient exploratory resampling methos. }  we need to take clustering into consideration.  \cite{Varoquaux12Clustering} proposed to first randomly rescaled the data matrix $\mathbf{X}$, then run the clustering during each iteration of stability selection, and then solve the corresponding problem (\ref{eq:everaged}). . Thus, the clustering will takes too much time, even prohibitive, especially for large scale problems. %If the stability selection repeats run in order to get different clusering
%\comm{  Challenge: Need to solve the sparse learning model NM times, where N: bootstrap samples, and M: number of parameter values.
% Efficient algorithms for sparse learning are needed. } Clustering is costing a lot of time.
%We propose to use the block resampling, which is efficient and also considers the grouping and strong correlation of voxels, to replace the randomize rescaled %clustering.  The blocked resampling is also even more efficient to shape the bounds of the discriminative regions. % The clustering helps to reduce the %correlation and reduce the searching space.

%Notice that in our paper, in order to further improve the stability of the result, we prefer to solve the model (\ref{eq:everaged}) based on averaging, rather than the more sophisticated model (\ref{eq:structuralL1}). Due to the very limited number of training samples, model  (\ref{eq:everaged}) is more reasonable in our cases, where we only need support information, i.e. the position of the involved voxels. %In addition, we further propose to use  the multi-stage to further reduce dimension of the voxel space to enhance the stability of results and the efficiency of our algorithm. %the number of resamplings and the expense of each clustering


%Therefore, due to the structured clusters and block sampling, it is possible for us to select more number of voxels than the number of control and patients. Here we are %caring about the number of clustered regions.


%\comm{}

%\item \comm{Why we can remove certain unrelated parts, to improve the accuracy of patterns? }
%\item Region average can also be implemented via PCA \cite{Margulies10fMRI}

%\item \comm{new resampling} spatial constrained resampling--named clustered resamplings
%\end{itemize}


%\subsection{Various Clustering Algorithms -- Hierarchical Clustering}
%\subsection{Various $\ell_1$ methods}

%\comm{Structural Sparsity}

%\comm{For now, only computational efficiency advantages. The theoretical verification is based on the previous study. about stability selection. The limitation on the ``laziness" of classifiers. }
%The above idea share something in common with the well known Recursive Feature Elimination \cite{Guyon02SVMRFE}, which is very effective in the field of gene selection. It is expected to  eliminates gene redundancy automatically and yields better and more compact gene subsets.
%However, notice that feature elimination might be not necessarily performing well consistently. For example, using the synthetic data introduce in the following numerical experiment section, we have observed the SVM-RFE performs not well, i.e. missing a significant portion of discriminative regions. A simple explanation is  \comm{if we want to make use of multi-stage, we need to }

%In this paper, we emphasize the tradeoff between the statistical inference quality and the computational efficiency, which is a significant advantage of our algorithm and makes it % is to find out the comprehensive involved brain regions in a computationally efficient way, especially
%suitable for very high dimensional voxel space.  The time costing mainly on the very small size sparse logistic regressions and very few number of clusterings.
%As demonstrated in numerical experiments,  it achieves much better computational efficiency compared with the common stability selection.  The procedure of our algorithm is summarized in following table. %  }





%\subsection{Multi-center Data}

%\subsection{Brief Theoretical Analysis and Support}
\section{Numerical Experiments} \label{Sec:Experiments}
In this paper, we will demonstrate
%compare our algorithm with the classical univariate voxel selection method, the recently popular
% state-of-the art multi-voxel pattern recognition methods including  ttest, $\ell_2$-SVM, $\ell_2$ Logistic Regression, and $\ell_1$-SVM, $\ell_1$ Logistic Regression, Elastic Net, and randomized  $\ell_1$ logistic regression. We have showed that while our method have achieved better control of false positives, it also has smaller  false negatives than other alternatives in \cite{Wang14RSS}. \comm{Here we focus more on the
the stability and reliability of the RSS based support identification method. This point will be empirically verified  based on the multi-center data, which is of data variation among different centers. In this paper, the data comes from
% \comm{Multi-modal} %We do not compare with recently randomized multivoxel method proposed in \cite{Varoquaux12Clustering}, because for our large scale problems, their algorithm is taking too much computing time due to too many times of run clustering algorithms and is not  feasible any more.  In fact,  our algorithm already achieves almost the best support identification performance one algorithm can achieve in our numerical experiments, in a computationally very efficient way.
%
% Notice that while there have existed similar comparison \cite{Li09VoxelSelection}, our comparison is among the first comprehensive ones, to our best knowledge, to  %\comm{Why not other popular methods?}
%
%One feature of our numerical experiments is as follows.  We use our carefully design synthetic data to
%demonstrate the limitation of the univariate methods and common multivariate methods in terms of support identification, the important role of randomization to overcome the ``laziness" of common classifiers, and the unique  role of the block subsampling in terms of reducing both false positive and false negative rate. %This point is often omitted in the previous research. Moreover, we compared with the common stability selection to demonstrate the better computational efficiency. %a very clear way.
%The true support is unknown for real data, easier to explain
%scores obtained on the best voxels outperforms alternative
%supervised methods suggesting that our approach
%is not only better in terms of support identification but
%can also significantly improve prediction on such data.
%
%There are differences between the $\ell_1$-norm logistic regression model
%(\ref{eq:L1}) and $\ell_1$-norm SVM. For instance, when $\ell_1$-norm SVM and
%model (2) are used for the same dataset, there are more variables
%and constraints for $\ell_1$-norm SVM than for (\ref{eq:L1}). This implies
%a heavier computational burden is required for $\ell_1$-norm SVM.
%Our algorithmic framework is also workable for $\ell_1$-norm SVM. This is an important feature of our new algorithmic framework,
%-expandability. The idea can be applied many other sparsity based classifiers to better identifier the discriminative brain regions.
%
%\subsection{data}
  the Autism Brain Imaging Data Exchange (ABIDE)- a grassroots consortium aggregating and openly sharing $1112$
existing resting-state functional magnetic resonance imaging (R-fMRI) data sets with corresponding structural MRI (sMRI) and phenotypic
information from 539 individuals with ASDs and 573 age-matched typical controls (TCs; 7-64 years) (http://fcon\_1000.projects.
nitrc.org/indi/abide/).  Then we follow the sample selection principle in \cite{Martino13MP} %, where imaging analyses were limited to: (1) male subjects, as
%they represented 90\% of the aggregate sample; (2) sites with IQ estimated
%for $>$ 75\% per diagnostic group; (3) individuals with FIQ within 2 s.d. of the
%overall ABIDE sample mean ($108\pm 15$); (4) individuals with mean framewise
%displacement (FD) $<1$ s.d. above the sample mean ($0.23\pm 0.27$ mm); (5) data with anatomical images providing near full
%brain coverage and successful registration; and (6) sites with at least 10
%participants per group after the above exclusions. This way, we have
get the final data for $763$
individuals (ASDs=360; TCs=403) from $14$ centers. %These multi-center data helps verify the reliability of support identification of our algorithm. %, which can then be applied to other related fields such as revealing the activated brain regions in terms of external multimedia stimulus.

{We compare with two sample t-test, $\ell_2$ logistic regression, $\ell_2$ SVM, TV-L1 \cite{Gramfort:2013:IPR:2552484.2552506}, $\ell_1$ logistic regression and randomized $\ell_1$ logistic regression (stability selection based on $\ell_1$ logistic regression). Two sample t-test and $\ell_2$-SVM are implemented as internal functions of MATLAB. The python code of TV-L1 is provided by Prof. Alexandre Gramfort of Telecom ParisTech and it is  under integration in the Nilearn package. Most of the other  algorithms have been implemented in LIBLINEAR \cite{REF08a} or SLEP (Sparse Learning with Efficient Projections) software \cite{Liu09:SLEP:manual}.} % except the two sample t-test and $\ell_2$-SVM, which have been implemented as  . }%Ones can check their detailed formulas in the manus of these softwares. %\comm{This part will be rewritten in more details}  % These methods

\subsection{Image preprocessing}
R-fMRI scans were preprocessed with ��SPM��(\url{www.fil.ion.ucl.ac.uk/spm/}). Image preprocessing steps
included slice-timing and motion correction, smoothing, correted R-fMRI
measures were normalized to Montreal Neurological Institute (MNI) 152
stereotactic space (3$mm^3$ isotropic) with linear and non-linear registrations.Linear detrend and temporal filter (0.01-0.08Hz) were then applied on the normalized images. For each participant, we generate the following voxel-wise regional metrics: Amplitude of Low Frequency Fluctuations (ALFF) \cite{YuFeng2007ALFF}.  ALFF is defined as the mean square root of the power spectrum
density over the low frequency band (usually 0.01,0.08 Hz). After
an individual ALFF map is obtained, %a standardization procedure
%can be performed that the individual
%ALFF map is divided by the mean ALFF value of voxels within a
%specific mask (whole brain mask or user-defined mask).
%While ALFF reflects the intensity of regional
%spontaneous brain activity, it is sensitive to the
%physiological noise found in some brain regions, such as in the
%cistern areas. An improved ALFF approach, namely,
%zALFF, is implemented in REST.
zALFF is obtained by performing  a standard z-transformation  on the voxels of individual ALFF map within a specific mask (gray matter mask). % to better suppress the noise.





For structural MRI (3DT1, here), we perform Voxel-based Morphometry (VBM), which involves a voxel-wise comparison of the local concentration of gray matter between two groups of subjects \cite{Ashburner20008VBM}.  All the skull-stripped and reoriented images were spatially normalized to the Montreal Neurological Institute (MNI) space by minimizing the residual sum of squared differences between structural MRI and the ICBM 152 template image. The data were then resampled to $3\times 3\times 3mm^3$. All these images were segmented into GM (grey matter), WM (white matter) and CSF (cerebrospinal fluid) using the unified segmentation algorithm with incorporated bias correction. GM images were then smoothed with an 8-mm smoothing kernel.  In this paper, we use the preprocessed VBM data as a main feature to test our algorithm.



%\subsection{Regional measures}
%For sMRI, voxel-based morphometry (VBM) involves a voxel-wise comparison of the local concentration of gray matter between two groups of subjects. The procedure is relatively straightforward and involves spatially normalizing high-resolution images from all the subjects in the study into the same stereotactic space. This is followed by segmenting the gray matter from the spatially normalized images and smoothing the gray-matter segments..

%30,31 a
%measure of the connectome graph indexing the number of direct
%connections for a given node.




\subsection{Settings of Algorithms}
 For our algorithm, %the blocked subsampling is proposed to be applied to the common multivariate classifiers such as  the sparse logistic regression (\ref{eq:everaged}).
 we use the  SLEP \cite{Liu09:SLEP:manual} software to efficiently solve the model (\ref{eq:everaged}), which in fact  a common $\ell_1$ model.
For the selection of regularization parameters of the involved multivariate methods except our method,  cross validation is used.
{ We set the subsampling rates are $0.5$ and $0.1$ for samples and voxels, respectivley. The resampling times is  $200$. These settings are following the common stability selection default settings \cite{Meinshausen10StabilitySelection}. }

% { We set  $N=1, K=200$ and correspondingly the  total resampling times is $N\times K=200$. That is to say, we are using a one-stage version. }

 {
  The internal k-means function of MATLAB is used for clustering. The default settings are used. For example, Squared Euclidean distance is adopted. One key parameter is the number of clusters, which is need to be prescribed.  The number of parameters are set according to the number of samples. For this multi-center data, each center has around $55$ samples. In general, we heuristically set the number of clusters between twice and $5$ times of the number of samples. We have tried $3$ different  numbers of clustering such as 116, 160 and 200, respectively,  in order to test how sensitive our final results are with this number. } %their effects on sensitive for .
 % for model selection,

{The block size is another key parameter of our algorithm.  It can be set according to the probably available prior knowledge of the size of the discriminative regions. In case of no such prior knowledge, we can give it a moderate value no larger than $10\times 10\times 10$. In the following subsection, we try different size of block from $3\times 3 \times 3$ to  $7 \times 7 \times 7$, which are considered to be in a reasonable range based on our experimental experiences. We can see that the final voxel selection result is not sensitive to the block size in this experiment. }%, we set it to be $6 \times 6 \times 6$. How sensitive it is? }, but still achieves an impressive performance.
%Previous fMRI study of large-scale brain network of board game experts identified significantly active brain regions in domain-related task. The case-control analysis (game condition vs. blank board condition) results indicated Task-Evoked Activation mainly in dorsolateral prefrontal cortex and parietal cortex of cognitive network \cite{Duan12chess}.
%Now we performed multivariate pattern analysis methods to identify discriminative brain functional activation regions in board game experts and compared the feature selection power of them.



{In order to control the false positives, we need to set a threshold value to filter out uninformative voxels, i.e.,   voxels whose corresponding weights have smaller magnitude than this threshold will be considered as noisy voxels. %The rest uninformative voxels will be not plotted.
The setting of a threshold value is quite difficult in general and may adopt different schemes in different situations. %A popular  method is to use the cross-validation based on prediction accuracy. The threshold value is chosen corresponding to the highest prediction accuracy. However, it is not suitable for multi-center data. }
%
%For fair comparison of different feature selection methods, we try different combinations of setting number of selected voxels for all the involved methods.
For two sample t-test method,  the number is generally determined by setting the p-value $ <0.05$ as significant level. %in  multiple hypothesis testing to correct for multiple comparisons, as commonly done. ttest�ǵ������ķ�����fdrУ��0.05����û��ͨ���ģ����ڲ�ȡ����ttest ��pֵС��0.05�����أ����о���ttest �� our method ��ȡ�˲�ͬ��cluster ��116 160 200���ķ����������l2logistic l2svm randL1s������ cluster�ĸ�����%and the same number features in different methods to summarize the results.
For  the multivariate methods, %we need to set the threshold value  for controlling the false positives, where the voxels of  corresponding weight is larger than the threshold value are selected.
%
%As mentioned before,
while the cross-validation
method is widely used, it usually works well for the cases where there is a lot of training data available
and the prediction accuracy is the main concern. In addition, it often causes large false negative rate though a small false positive rate can be often achieved. % Here we would like to have a small false negative rate for the mutual validity of selected features among different centers.
} In this paper, we use a very simple but flexible way to set the threshold value.  It is a probability based method suggested in \cite{Li09SparseVoxelSelection}.  We repeat the description of the probability method in \cite{Li09SparseVoxelSelection} as follows.  %The threshold value is determined
%by the cross-validation method and the probability method
%are generally different. Here we adopt method suggested in \citep{Li09SparseVoxelSelection}.
%
%
Specifically, considering that the entries of $\mathbf{w}$ are sparse,
we assume that the probability distribution of the entries of
$\mathbf{w}$ is Laplacian. Using all entries of $\mathbf{w}$ as samples, we estimate
the mean, the variance, and the inverse cumulative distribution
function F-1 of this Laplacian distribution. We then
define $R = \{i: |\mathbf{w}_i| > \theta, i = 1, �� �� ��,p\}$, where $\theta$ is chosen as
F-1 ($p_0$), $p_0$ is a given probability (e.g., $0.975$ or more restrictive $0.99$ in this paper). %\comm{Why this value}.
{  Unlike the single center data analysis, the threshold is not necessarily very  restrict because  the operation of intersection of the selected features between different sits helps remove false positives.  } Moreover,  we can gradually reduce the $p_0$ from a large value (for example, $0.99$) and check the corresponding selected voxels, which are the intersection of the selected voxels of each center. We have observed that when $p_0$ is reduced to a certain value (for example, 0.95), the selected voxels would  stay unchanged. In such cases, we have obtained a reasonable threshold. We will show this heuristic rule in the following description of the numerical results.
%multi-center can help verify the true positives and remove the false estimations. }




 \subsection{Evaluation Criteria}
% The main feature of our algorithm is the improved stability and sensitivity.
We would like to  demonstrate that our method can achieve  a consistent feature selection result across different centers, due to  our  incorporation of  ideas of stability selection and structural sparsity.
%
%{Threshold is not necessarily strict, for multi-center data. }
%
 Since there is no ground truth for  evaluation, the confidence level of a selected relevant voxel depends on the number of centers where it is commonly selected. Since the feature selection algorithm is performed on each center  individually, a voxel appearing at the voxel selection results of many centers is more likely to be true informative voxel. Therefore, we will show and compare the consistent selected voxels among different centers of all these involved algorithms.

{ In addition, in order to further demonstrate that  these intersected regions discovered  by our algorithm are  stable, we  perform a false positive estimation scheme based on the permutation test and cross validation, which is proposed in \cite{Rondina14}. It  aims to calculate the ratio of false positives among all the finally selected voxels. % , which is based on a permutation test and cross validation.

Finally, we also present the prediction power of the selected voxels. The details are in the following subsection \ref{subsec:analysis}.} %\comm{Notice that }

\subsection{Clustering Results and Feature Selection Results}
\begin{figure*}[!h]
    \centering
    \includegraphics[width=0.7\textwidth,height=0.18\textwidth]{SDSU_fMRI_kmeans.eps}\vspace{0.1cm}
    \includegraphics[width=0.7\textwidth,height=0.18\textwidth]{SDSU_fMRI_results.eps}
    \vspace{-0.1cm}\caption{The upper row is the clustering result by K-means with the number of clusters being $200$, for the fMRI data of SDSU center. The lower row is the brain map of the selected voxels on the same center, with the threshold parameter $p_0=0.975$.
   }\label{Fig:Clustering_fMRI}
\end{figure*}

\begin{figure*}[!h]
    \centering
    \includegraphics[width=0.7\textwidth,height=0.18\textwidth]{NYU_sMRI_kmeans.eps}\vspace{0.1cm}
    \includegraphics[width=0.7\textwidth,height=0.18\textwidth]{NYU_sMRI_results.eps}
    \vspace{-0.1cm}\caption{The upper row is the clustering Results by Kmeans with the number of clusters being $200$, for the sMRI data of NYU center. The lower row is the brain map of the  selected voxels on the same center, with the threshold parameter $p_0=0.975$.
   }\label{Fig:Clustering_sMRI}
\end{figure*}



{As mentioned before, the proposed RSS challenged the limitations of stability selection. The original stability selection  fails to explicitly make use of the assumption that these targeted voxels are often spatially contiguous and results in distributed brain regions. We used K-means to cluster the voxels into groups to model the spatial constraints. We need to emphasize that  the clusters resulted from K-means can not accurately model the underline homogeneity of the brain regions. However, due to the adoption of subsamplings, we can reduce the bias of the final feature selection results caused by K-means-based clustering.  Figure \ref{Fig:Clustering_fMRI} is the K-means clustering result and corresponding feature selection result on fMRI data of the SDSU center.
We can see the results of K-means clustering are reasonable and roughly match the prior  knowledge of brain segmentation and integration to some degree. In addition, while the clustering information helps reduce the variance, the final selected brain regions do not necessarily exactly follow the grouping based on the K-means clustering. This way, the possible bias caused by the clustering is expected to be reduced. We can also observe similar phenomenon from  sMRI data, as showed in  Figure \ref{Fig:Clustering_sMRI}.  }

{  We have also tested $3$ different  numbers of clusters such as $116$, $160$ and $200$. Figure \ref{Fig:Our1-clusterignNumberSMRI} and Figure \ref{Fig:Our1-clusterignNumberFMRI} are the results of our algorithm applied to sMRI and R-fMRI data, respectively. We have observed that the final support identification results are not very  sensitive to the number of clusters. The revealed brain regions are the same, though they might be of slightly different size for different settings. }
\begin{figure*}[!h]
    \centering
    \includegraphics[width=0.8\textwidth,height=0.85\textwidth]{T1_overlapped_new.eps}
    \vspace{-0.1cm}\caption{Selected overlapped voxels of our method shared by all 14 centers for the sMRI, i.e., 3DT1 data. We tested different settings of number of clusters (clusters no. ) and different $p_0$ (prob.level) to demonstrate the robustness of our algorithm.
   }\label{Fig:Our1-clusterignNumberSMRI}
\end{figure*}


\begin{figure*}[!h]
    \centering
    \includegraphics[width=0.8\textwidth,height=0.85\textwidth]{zalff_overlapped_new.eps}
    \vspace{-0.1cm}\caption{Selected overlapped voxels of our method shared by all 14 centers for the fMRI data. We tested different settings of number of clusters (clusters no. ) and different $p_0$ (prob.level) to demonstrate the robustness of our algorithm.
   }\label{Fig:Our1-clusterignNumberFMRI}
\end{figure*}

%\comm{How is the clustering results are interpretable}
\subsection{{ Feature Selection Results Corresponding to Different Block Sizes}}
{Now we are checking the robustness of our algorithm in terms of different block sizes when performing block subsampling. While we do not have the formula for  the optimal block size,  we show that the final feature selection result is not sensitive to the block size. We tested $5$ different block sizes, i.e. $3\times 3 \times 3=27$, $4 \times 4 \times 4=64$, $5  \times 5 \times =125$, $6 \times 6 \times 6=216$ and $7 \times 7 \times 7=343$. }  Figure \ref{Fig:Our2-blocksize-fMRI}  is the feature selection result when our algorithm is applied to the multi-center fMRI data set. We can see that  the results corresponding to different block sizes only have a small difference in term of the size of detected areas. Moreover, they all indicate the same  brain locations.  The similar conclusion can also be drawn for sMRI data, as suggested by Figure \ref{Fig:Our-blocksize-sMRI}.

\begin{figure*}[!h]
    \centering
    \includegraphics[width=0.8\textwidth,height=0.75\textwidth]{zalff_overlapped_different_blocksize_0975}\vspace{0.1cm}
    \vspace{-0.1cm}\caption{Feature selection results of our algorithm corresponding to different block sizes for fMRI data.
   }\label{Fig:Our2-blocksize-fMRI}
\end{figure*}

\begin{figure*}[!h]
    \centering
    \includegraphics[width=0.8\textwidth,height=0.75\textwidth]{T1_overlapped_different_blocksize_099}\vspace{0.1cm}
    \vspace{-0.1cm}\caption{Feature selection results of our algorithm corresponding to different block sizes for sMRI data
   }\label{Fig:Our-blocksize-sMRI}
\end{figure*}

\subsection{Comparison of Consistency of Results of Different Algorithms}

%\subsection{\comm{Comparison of Consistency of Results of Different Algorithms}}
%{\footnote
\begin{table}[!h]\tiny
\begin{tabular}{|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l}
  \hline
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
% &  SO-SA &DYCORS&	LMSRBF	& DDS	&SSKm	&ESGRBF	&NOMADm	&EGO\\\hline
%Rastrigin5&  -4.99&	-4.28	&-4.88	&-4.96	&-4.38	&-4.24	&-4.98\\\hline
Algorithms&ttest&	\multicolumn{2}{c|}{L2logistic}	&\multicolumn{2}{c|}{L2SVM}	&\multicolumn{2}{c|}{L1Logistic}		&\multicolumn{1}{c|}{TV-L1}&\multicolumn{2}{c|}{Randomized L1}&\multicolumn{6}{c|}{Our method}\\\hline
No. Clusters& &	\multicolumn{2}{c|}{}	&\multicolumn{2}{c|}{}	&\multicolumn{2}{c|}{}&\multicolumn{1}{c|}{}	&\multicolumn{2}{c|}{}&\multicolumn{2}{c|}{116}&\multicolumn{2}{c|}{160}&\multicolumn{2}{c|}{200}\\\hline
prob.level &	& 0.99    &0.975	& 0.99&0.975	& 0.99&0.975& & 0.99&0.975&0.99&0.975& 0.99&0.975& 0.99&0.975\\\hline												KKI	&	    1303&		10873&	16438&		10454&	16337&	12	&17&  8290&	8491&	14274&			8583&	8583&	8459&	8459	&8146	&8146\\\hline
Leuven&		1590&		10440&	15426&		10788&	16517&	51	&51	&	6963&8177&	13266&			7699	&7699&	7753	&7753	&8487&	8487\\\hline
MaxMun	&	1983&		10456&	15511&		11537&	17568&	63&	64	&	7492& 8250&	13145&			8105&	8105	&7975&	7975	&8395&	8395\\\hline
NYU	&	    2160&		10856&	16327&		10687&	16395&	422&455&	7405& 8539&	14572&			8109&	8109&	8108&	8108	&8334&	8334\\\hline
OHSU&		1643&		9508&	14368&		11158&	17148&	12&	12	&	6737&8009&	12853&			9007&	9007&	9105&	9105&	9693&	9693\\\hline
Olin&		958	&	    10259&	15196	&	11377&	17171&	1&	1	&	6332&8260&	13343&			9013&	9013&	8932&	8932&	8777&	8777\\\hline
Pitt&		954	&	    10764&	16084	&	11398&	17298&	8&	8	&	5341&8562&	14355&			8593&	8593&	8632&	8632&	8678&	8678\\\hline
SDSU&		1915&		10648&	15590	&	11970&	17985&	6&	6	&	7176&8628&	13362&			9766&	9766&	10477&	10477&	9753&	9753\\\hline
Standford&	2457&		10970&	16309	&	11061&	16869&	50&	50	&	10642&8655&	13981&			9063	&9063&	9033	&9033&	8636	&8636\\\hline
Trinity		&2423&		10246&	15355	&	11351&	17556&	9	&9	&	9370&8270&	12911&			8424	&8424	&8297&	8297&	8695&	8695\\\hline
UCLA	&	23643&		10659&	15921	&	10952&	16746&	178	&184&	7934&8324&	13729&			7547&	7547&	7807&	7807&	8092&	8092\\\hline
UM	&	    2587	&	10793&	16237	&	11571&	17607&	387	&413&	7043&8628&	14212&			8230&	8230&	8621&	8621&	8855&	8855\\\hline
USM	&	    8298&		10523&	15402	&	11577&	17181&	323	&330&	7714&8393&	13738&			7896&	7896&	7952&	7952&	7767&	7767\\\hline
Yale&		3801&		10462&	15617	&	11201&	17275&	2	&2	&	5749&8312&	12916&			9143&	9143&	8862&	8862&	8592&	8592\\\hline
overlapped(14)&		0&		0&	   2	&	0&	0&	0&	0&		0&0&	1&			1436&	1436	&1586&	1586&	1432	&1432\\\hline
overlapped(11)	&	0&		75&	531	&	0&	135&	0&	0&		0&177&	1664&			4557	&4557&	4674&	4674	&4838&	4838\\\hline
overlapped(9)	&	0&		703&	3430&		152&	2133&	0&	0&	0&	1101&	5924	&		6420&	6420&	6526&	6526&	6643	&6643\\\hline
%overlapped(5)	&	70&		12894&	25234&		11433&	28304&	8&	8&	1015&	13823&	25714	&	10736	&	10736&	10552&	10552&	10756	&10756\\\hline

\end{tabular}
\caption{Summary of Number of Selected Voxels for the multi-center sMRI data}\label{tab:summary2}
\end{table}

%{\footnote
\begin{table*}[!h]\tiny
\begin{tabular}{|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l}
  \hline
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
% &  SO-SA &DYCORS&	LMSRBF	& DDS	&SSKm	&ESGRBF	&NOMADm	&EGO\\\hline
%Rastrigin5&  -4.99&	-4.28	&-4.88	&-4.96	&-4.38	&-4.24	&-4.98\\\hline
Algs&ttest&	\multicolumn{2}{c|}{L2logistic}	&\multicolumn{2}{c|}{L2SVM}	&\multicolumn{2}{c|}{L1Logistic}	& \multicolumn{1}{c|}{TVL1} &\multicolumn{2}{c|}{Randomized L1}&\multicolumn{9}{c|}{Our method}\\\hline
No. & &	\multicolumn{2}{c|}{}	&\multicolumn{2}{c|}{}	&\multicolumn{2}{c|}{}	&\multicolumn{1}{c|}{}&\multicolumn{2}{c|}{}&\multicolumn{3}{c|}{116}&\multicolumn{3}{c|}{160}&\multicolumn{3}{c|}{200}\\\hline
prob.level &	& 0.99    &0.975	& 0.99&0.975	& 0.99&0.975& & 0.99&0.975&0.99&0.975&0.96& 0.99&0.975& 0.96 & 0.99&0.975&0.96\\\hline												
KKI		&4347	&	10565&	15683&		11096&	16510&	249	&258  &5699  &	8380&	13286	&		9174&	9174& 9174&	7448	&7448&	7448 &7652	&7652&7652\\\hline		
Leuven	&2485	&	9884&	15090&		11709&	17578&	261&273	  &5868  &	8079&	13027	&		8293&	15210&	15210& 11964&	14248	&14248	&12720&14998&14998\\\hline		
MaxMun	&2819	&	9746&	14667&		11523&	17317&	277	&286 &5994&		8022&	12723	&		8872&	8872& 8872&	10148&	10148&	10148	&9964&9964&9964\\\hline		
NYU		&8339	&	9928&	15051&		11286&	16750&	653	&696	&5604&	8153&	13480	&		8290&	8290&8290&	8055&	8055&	8055	&8181&8181&8181\\\hline		
OHSU	&1785	&	9214&	13853&		11985&	17808&	239	&248	&5938&	7837&	12300	&		10831&	10831&	10831& 12604&	12604&	12604	&8077&12040&12040\\\hline		
Olin	&1532	&	10534&	15656&		11639&	17377&	189	&193	&5252&	8252&	13371	&		9618&	9618&	9618& 9515&	11491&	11491&	7981&11364&11364\\\hline		
Pitt	&1988	&	9966	&15057&		11847&	17811&	296	&321	&5863&	7999&	12901	&		8522&	8522& 8522&	7930&	8229&	8229&	8726&8726&8726\\\hline		
SDSU	&3539	&	10185&	15189	&	12117&	18201&	222	&229	&5938&	8635&	13629	&		11461&	11461& 11461&	10435&	10435&	10435&	9966&10412&10412\\\hline		
Standford	&	2486&10806&	16059	&	11772&	17467&	211	&217	&5853&	8529&	13836	&		9616&	9616& 9616&	9400&	9400&	9400&	11640&11640&11640\\\hline		
Trinity		&2742	&9874&	14815	&	12164&	18045&	349	&369	&5598&	8104&	12848	&		8385&	8385& 8385&	8964	&8964&	8964	&8447&8447&8447\\\hline		
UCLA		&1166	&9703&	14786	&	11818&	17501&	486	&500	&5516&	8046&	13191	&		8507&	10002&10002&	8895&	9896	&9896&	7920&9854&9854\\\hline		
UM		&3545	&10645& 	16123	&	11556&	17277&	508	&537	&5706&	8423&	13890	&		9417&	12500&12500&	10066&	12204&	12204&	8679&11344&11344\\\hline		
USM		&5238	&	9618&	14456	&	12003&	17765&	488	&511	&6006&	8183&	12954	&		7359&	7359&7359&	8892&	8892&	8892&	8786&8786&8786\\\hline		
Yale	&	2907&9112&  	13684	&	11314&	16652&	175	&180	&5874&	7609&	12026	&		11033&	11033&11033&	10598&	10598&	10598&	11630&11630&11630\\\hline		
OL(14)	&	0&	0&	2	&	0      & 	   1&	0	&0	&0	&1&	24	&	645 &	721&	721	& 688	&771&771	& 801&857&857\\\hline		
OL(11)	&	0&	104&	583&		4&	153  &	0	&0	&	0&158& 	1129& 2708&	3118&	3118&	2946&	3224&3224&	3131&3596&3596\\\hline		
OL(9)	&	0&	872&	3344&		153&2298	&0	&0	&	0&892&	3667& 4558&	4871&	4871&	4896&	5180&5180	&4968&5427&5427\\\hline	
%OL(5)	&	88&	11216&	22395&		13049&30908	&1	&1	&	1012&7770&	16980&10075&	11080&	11080&	10752&	11546&	11546&10181&11828&11828\\\hline	
\end{tabular}
\caption{Summary of Number of Selected Voxels for the multi-center fMRI data}\label{tab:summary1}
\end{table*}
%}



  All the univariate and multivariate pattern feature selection methods can identified certain number of possible disease-related voxels for each center. The number has been listed in Tables  \ref{tab:summary2} and \ref{tab:summary1}, which corresponds to the results on the sMRI and fMRI, respectively.  Besides TV-L1, the multivariate methods mostly used  different $\theta$ values (denoted as prob.level), i.e. $0.99$, $0.975$ and even $0.96$ when applying the probability method to determine the selected voxels. For TV-L1, we have the third-party python software and we use its internal default way to set the threshold value, which is partially based on cross-validation.  The overlapped($S$), represents the number of the overlap of selected voxels across different centers, where $S$ means that at least $S$ centers share these selected voxels. As the $S$ increases, the number of overlapped voxels might decrease as expected.

  The  variation of MRI data between different centers brings both challenges and benefits.
    On one hand, this kind of variation makes many existing state-of-the-art algorithms likely to select quite distinct voxels for different centers,   because they tend to select a
minimum subset of features to construct a classifier of the best predictive accuracy and  often ignore
``stability" of feature selection \cite{Yu08denseFeatureGroups,He2010215,buhlmann2011controlling}. Therefore, the results across different centers could be quite different and the overlapped parts are quite few. From  Tables  \ref{tab:summary2} and \ref{tab:summary1} , we can see that the overlapped voxels corresponding to the two sample t-test, L2-logistic, L2-SVM, TV-L1, and L1-Logistic are almost none when $S$ is larger than $11$. Notice that while  two sample t-test is widely believed to be a stable method for single center data analysis, it is not a good choice for multi-center data. TV-L1, as a recently popular feature selection method for the single center data analysis, does not work well for the multi-centere data analysis, partially due to the lack of stability scheme.   As expected, randomized $\ell_1$ logistic regression, as a method of stability selection, behavior  better than other alternatives. However, it is still much worse than our method, especially when $S$ is large, for example, $11$ or $14$, partially due to the failure of explicitly making use of the prior structural information of the discriminative voxels. For better comparison, We also plotted the selected voxels of different algorithms corresponding to $S=5, 9, 11, 14$, in Figures  \ref{Fig:Compare9_T1}, \ref{Fig:Compare11_T1} and \ref{Fig:Compare14_T1}  for the sMRI data, and  in Figures  \ref{Fig:Compare9_fMRI}, \ref{Fig:Compare11_fMRI} and \ref{Fig:Compare14_fMRI} for the fMRI data, respectively. %Due to the space limitation, we did not plot the results of fMRI, though it shares almost the same conclusion.
We can see that our method can always achieve a consistent voxel selection result across different centers while the other state-of-the-arts mostly fail. %, though they could work  for a single center data analysis.
%
  %prefrontal and parietal lobe regions.

  On the other hand, the variation of data between different centers can help us to verify the reliability of the selected discriminative voxels \cite{Friedman08retest}.   If these selected voxels are all shared by these involved centers, a better reliability of them is expected, even as the potential biomarker in the possible clinical diagnosis. Despite the computationally challenge, our algorithm incorporating the ideas of stability selection and structural sparsity successfully revealed significant number of the shared selected voxels.  %Figures \ref{Fig:Our1} and \ref{Fig:Our2} correspond to the results of sMRI and fMRI, we will explain their reasonability and meanings from the disease point of view in the next section.

{Moreover, we would like to point out that  there exists a practical simple way to set the threshold value for feature selection for our algorithm. We have empirically observed that the sorted absolute values of the weights (from large to small) first gradually change then suddenly decreases dramatically and produce a gap. Usually we can set the threshold value at the location of the gap, and it is expected to achieve a good control of false positives. A similar  phenomenon has been observed and exploited in our previous related work about sparse signal recovery \cite{wang2010sparse}. Specifically, let us look at the Table \ref{tab:summary2}. For the two threshold values corresponding to probability levels $p_0=0.99$ and $p_0=0.975$, respectively,  the same discriminative voxels are obtained for our algorithm as each of the center. That is to say, the sorted weights of voxels indeed have a gap. Using this threshold value, we obtained an interpretable discriminative brain regions and we will give further analysis of these regions in the following subsection.  Table \ref{tab:summary1}  tells a similar story, except the gap is around the threshold values corresponding to the probability levels $p_0=0.975$ and $p_0=0.96$.  }


%fast and the value at the change point, surrounded by a red circle, is an appropriate value. This appearance of the inflexion point is not only because that our method is based on sparse regularization and therefore most of the components of the $\mathbf{w}$ are forced to be close to $0$, but also because that our method can achieve very high voxel selection accuracy for our synthetic data as showed in the last subfigure. Notice that the original stability selection (Randomized L1 logistic regression) does not have such an obvious threshold point and its has  worse specificity and sensitivity than our method as showed in the last subfigure. Using the threshold value suggested by the changing point, we achieve a good support recovery as showed in the fifth subfigure. The same number of top weighted voxels of the randomized $\ell_1$ logistic regression are plotted in the fourth subfigure, which mostly misses the first discriminative block and introduces many false positives. Notice that our ``constrained block subsampling" does not bring blurring to the final result, as also showed in the last subfigure.

%{The intersections between different sites helps remove the false positives. }

%\begin{itemize}
%\item Our method can obtain  bigger or many more discriminative regions.
%\item These regions are explainable.
%\item \comm{The detailed results will be given soon. }
%\end{itemize}

\subsection{{Subsequent Analysis of Feature Selection Results}} \label{subsec:analysis}
{We have observed that the overlap of selected voxels   by our algorithm is significantly higher. We believe that most of the selected voxols are true positives because the false positive usually are unlikely to repeat on every center. } Furthermore, we would like to further estimate the ratio of  false positives from another point of view based on permutation test and cross validation, as suggested in \cite{Rondina14}. The results are summarized in Table \ref{tab:summary-sMRI-fMRI}.  For sMRI, the estimated false positives ratios are all smaller than 3\%, for the overlap  of all the $14$ centers, for all the settings using different clustering numbers (116, 160, 200) and different probability levels (0.99 and 0.975). For fMRI, the estimated false positives ratios are slightly larger, but still smaller  than 6\%.  Notice that these numbers are only rough estimations for reference, not necessarily accurately reflecting the ground truth.  In addition, we provide the classification performances when using those identified voxels as features to construct a classifier. Notice  that for many alternative methods, they could not achieve consistent results if $S$ is large and so we have to reduce $S$ to $5$ or even $2$. For our method, we use $S=11$. When performing the classification, we pool all samples from $14$ centers together and we have in total $763$ samples. We randomly pick $700$ as the training samples and the rest $63$ as the test samples and we repeat this procedure for $100$ times. We use the $\ell_2$ logistic regression as the classifier. The classification accuracy is summarized in Tables \ref{tab:prediction-sMRI} and \ref{tab:prediction-fMRI} for sMRI data and fMRI data, respectively. While the selected voxels by our algorithm can achieve the best classification performance, they are generally not very high. A possible reason is that  a simple logistic classifier can not deal with the relatively high data variation of these samples from different centers and we might consider more sophisticated classifiers in the future.   %In addition, we need to m


%\comm{Not finished, will be added soon}

\begin{table}[!h]\tiny
\begin{tabular}{|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l}
  \hline
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
% &  SO-SA &DYCORS&	LMSRBF	& DDS	&SSKm	&ESGRBF	&NOMADm	&EGO\\\hline
%Rastrigin5&  -4.99&	-4.28	&-4.88	&-4.96	&-4.38	&-4.24	&-4.98\\\hline
&\multicolumn{6}{c|}{Our method}\\\hline
No. Clusters& \multicolumn{2}{c|}{116}&\multicolumn{2}{c|}{160}&\multicolumn{2}{c|}{200}\\\hline
prob.level &0.99&0.975& 0.99&0.975& 0.99&0.975\\\hline												
No. Totally Selected&					1436&	1436	&1586&	1586&	1432	&1432\\\hline
No. Estimated False Positives	&	25	& 25&	37&	37	&23&	23\\\hline
\end{tabular}\begin{tabular}{|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l}
  \hline
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
% &  SO-SA &DYCORS&	LMSRBF	& DDS	&SSKm	&ESGRBF	&NOMADm	&EGO\\\hline
%Rastrigin5&  -4.99&	-4.28	&-4.88	&-4.96	&-4.38	&-4.24	&-4.98\\\hline
&\multicolumn{6}{c|}{Our method}\\\hline
No. Clusters& \multicolumn{2}{c|}{116}&\multicolumn{2}{c|}{160}&\multicolumn{2}{c|}{200}\\\hline
prob.level &0.99&0.975& 0.99&0.975& 0.99&0.975\\\hline												
No. Totally Selected&					645&	721	& 688&	771&	801	&857\\\hline
No. Estimated False Positives	&	24	& 33&	27&	39	&45&	51\\\hline
\end{tabular}\caption{Number of Selected Voxels (across 14 centers) and Estimated False Positives. The left part corresponds to the  multi-center sMRI data and
 the right part corresponds to the  multi-center fMRI data}\label{tab:summary-sMRI-fMRI}
\end{table}


\begin{table}[!h]\tiny
\begin{tabular}{|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l}
  \hline
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
% &  SO-SA &DYCORS&	LMSRBF	& DDS	&SSKm	&ESGRBF	&NOMADm	&EGO\\\hline
%Rastrigin5&  -4.99&	-4.28	&-4.88	&-4.96	&-4.38	&-4.24	&-4.98\\\hline
%\multicolumn{7}{c|}{sMRI}\\\hline
ttest&l2logistic&L2SVm&L1Logistic&TVL1&Randomized L1&Ours\\\hline
overlapped(5)&overlapped(11)&overlapped(11)& overlapped(2)& overlapped(5) &overlapped(11)& overlapped(11)\\\hline
0.6046&0.6054&0.6121&0.6012&0.6113&0.5941&0.6332\\\hline
\end{tabular}\caption{Prediction power of the selected voxels for sMRI data}\label{tab:prediction-sMRI}
\end{table}

\begin{table}[!h]\tiny
\begin{tabular}{|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l|l}
  \hline
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
% &  SO-SA &DYCORS&	LMSRBF	& DDS	&SSKm	&ESGRBF	&NOMADm	&EGO\\\hline
%Rastrigin5&  -4.99&	-4.28	&-4.88	&-4.96	&-4.38	&-4.24	&-4.98\\\hline
%\multicolumn{7}{c|}{fMRI}\\\hline
ttest&l2logistic&L2SVm&L1Logistic&TVL1&Randomized L1&Ours\\\hline
overlapped(5)&overlapped(11)&overlapped(11)& overlapped(2)& overlapped(5) &overlapped(11)& overlapped(11)\\\hline
0.6062&0.5908&0.6145&0.5944&0.5997&0.5801&0.6283\\\hline
\end{tabular}\caption{Prediction power of the selected voxels for fMRI data}\label{tab:prediction-fMRI}
\end{table}

  %
   %  {False positive estimate and prediction accuracy} We will also give a false estimate test to the final selected voxels and consider their prediction powers.

%\comm{In particular, we use the selected voxels of different algorithms to construct different classifiers, respectively.  We compared the different classification accuracy of these classifiers. We also construct a classifier without voxel selection, i.e. using all the voxels. }

\subsection{Discussion Based on Biomarkers from Multimodal Data}
We aim to show that these discovered biomarkers are interpretable and correspondingly support the reliability and stability of our algorithm. %and give some explanations for ASD.
According to Figures \ref{Fig:Our1-clusterignNumberSMRI} and \ref{Fig:Our1-clusterignNumberFMRI} , we found the grey matter abnormalities in precuneus, posterior cingulate cortex, insular, middle temporal gyrus-R, fusiform gyrus, hippocampus, parahippocampal gyrus may have the potential be the related biomarkers of sMRI data(VBM) for distinguishing between the autism and the healthy controls. And with respect to the fMRI (zalff), medial prefrontal, precuneus, temporal pole, superior and medial frontal gyrus be the revealed regions as the discriminations of ASD.
The potential structural biomarker regions are mainly concentrated in the DMN (default mode network), such as precuneus and posterior cingulate cortex(PCC) which are regarded relating to primary responsible for autobiographical memory and self-reference processing. For example, precuneus and  PCC are primarily responsible for autobiographical memory and the self-reference processing. Significant abnormality in insula was also revealed, it is a critical component of SN(salience network) and would be involved in switching between central-executive and DMN \cite{Sridharan26082008}. And the SN is thought to regulate dynamic changes in other networks and the damage to the structure of SN should disrupt the regulation of associated networks such as DMN \cite{Bonnelle20032012}. Besides the brain regions mentioned above, some other regions such as fusiform gyrus, hippocampus and middle temporal gyrus-R�� which are reported in previous studies \cite{Schumann14072004}. In our research, those overlapped regions we identified as discriminatied features, are found across all the 14 centers. Thus we can say that they may be the ASD-related stable biomarkers for the diagnoses and further analysis, according to the characters of our algorithm purposed and discussed before.
Similar to the structural data, the functional abnormalities of the ASD mainly lie in regions of medial prefrontal (mPFC), precuneus and temporal pole (above three are the components of the DMN) and anterior cingulate(ACC). mPFC is in charge of social cognition associated with self and others \cite{Gusnard27032001}.  In addition, the alterations of ACC, a part of the brain��s limbic system, may be relevant to the character of autism \cite{Bush2000215}. Accordingly, the consistent  overlapped features in ASD from multi-center data further confirmed the functional abnormalities in ASD \cite{Castelli02,Critchley00}.


  %However, besides each of the above methods has big false negative rate, meaning that it is only able to find out part of the comprehensive discriminative %regions. In addition, due to the $\ell_2$ regularization encouraging too much continuousness, not simplicity or sparsity, these results have significant %false positives.%}

\section{Conclusion and Future work} \label{Sec:conclusion}
%The success of our approach with very small sample
%sizes lies in exploiting the characteristics of the problem tackled: recovery of contiguous clusters from image
%measurements with local correlations. We reduced the number of clustering by the blocked sampling techniques and
%multistage refineness procedure.





% if have a single appendix:
%\appendix[Proof of the Zonklar Equations]
% or
%\appendix  % for no appendix heading
% do not use \section anymore after \appendix, only \section*
% is possibly needed

% use appendices with more than one appendix
% then use \section to start each appendix
% you must declare a \section before using any
% \subsection or using \label (\appendices by itself
% starts a section numbered zero.)
%
In this paper, we consider a commonly existent problem in mental disease diagnosis or congtive science based on brain imaging, i.e. locating the discriminative or activitaed brain regions correspoing to the specific decease or certain external stimulus including looking at multimedia materials. We revisit our recently proposed algorithm which incorporates the ideas of stability selection and structural sparsity, and aims to show its abilitity to reliably and stably find out these intrested brain regions  despite the inherent data variations or noise.
We consider to  demonstrate our point via the multi-center MRI data analysis, which is also a promising research direction by itself.  Specifically, our algorithm helps us to discover a set of consistent selected voxels across different centers and these voxels are quite interpretable and can be possible biomarkers for the correspoinding medical diagnosis. Considering that most of existing algorithms fail to overcome the difficulty brought by the variation between centers and therefore are not able to get a consistent result, our algorithm also seems to be promising for the multi-center data analysis from the algorithmic point of view. However, some therotical analysis of the performance of our algorithm will be required, considering most of current study are more on the empirical side.
%analysis problem. While  most of existing research has focused primarily on establishing results for prediction consistency, we focus more on the quality of feature selection, since reliable variable selection is important in many scientific areas. % where simple models are needed to provide insights into complex systems
%For the high dimensional fMRI data considered in this paper, %voxel selection is very important for both encoding and decoding of the measured data. %For decoding, it helps reduce the overfitting by dimensional reduction. As
%the encoding, our new method give a more comprehensive activity map compared with other existing state-of-the-art algorithms.
%Correspondingly,
%we propose a simple and computationally efficient approach for data-driven
%oxel selection, which is also called support identification potentially for the following biomarker extraction \cite{Orru2012SVM} in different literature. We propose to use the ``constrained block subsampling" as a variant of the classical stability selection, applied to the existing sparse multi-variate classifiers such as $\ell_1$ logistic regression, in case of fMRI data, which has the strong correlations and distributed multivariate discriminative patterns.
%, because the block subsampling corresponds to the multivariate pattern, which is reflected in certain number of neighboring voxels consisting of a discriminative brain region.
%Theoretically, we need to consider all possible combinations of groups, but this is almost impossible in practice.
 %For better computational efficiency and inference quality, we can take
  %  data driven clustering (as did in this paper) or prior anatomical information into consideration and this is where the ``constrained" comes from.  %The above exploration procedure is realized via the randomized subsampling in this paper and is able to accurately shape out all the discriminative regions in our cases.  In addition,
    %The   interpretability and stability of the resulted sparse and localized compact solutions are quite relevant to providing a credible diagnostic \cite{jenatton12multiscale}.
%
  %While our algorithm can be considered as a multivariate pattern recognition based method,
%It can also overcome the "laziness" of these methods in terms of reducing the false negatives.
 %In addition, the introduced block subsampling make it more efficient and effective than classical stability selection in terms of reduced false negative rate.
% The block subsampling corresponds to the multivariate pattern, which is reflected in certain number of neighboring voxels consisting of
%a discriminative brain region.
%Numerical experiments demonstrate that our method can achieve much better both specificity and sensitivity than other alternatives in terms of voxel selection accuracy rather than prediction or classification accuracy.
%Notice that while we take the sparse logistic regression as an example in our paper, the idea of randomized block subsampling can be applied to other sparse multivariate  pattern recognition methods such as the $\ell_1$-SVM. We evaluate the voxel selection quality of different methods via ROC in terms of
%
%In this paper, we are currently only focusing on the voxel selection based on the training data.
%
%Considering that the accurate  identification of discriminative voxels helps produce a classifier with better prediction accuracy on the test data in most cases \cite{DeMartino200844},
 %
 We will further test the performance of our algorithm through  more kinds of multi-center data and feature types
 %in terms of the classification accuracy on the test data when combined with common classifiers
 in the future work, epecially in the cognitive experiments including the study of brain responses to multimedia.
%{We have not yet had a good classifier for mutli-center data.}

%{bias variance tradeoff--using a more uniform template}
%\comm{Finally, the sample size used here was small, and although this was sufficient for these purposes, a larger sample size may be required for future studies, particularly if more complex tasks are administered that might evoke more subtle neuronal responses.}

 %to our algorithm can find out

%We emphasize the exploration aspect during the detection of  the discriminative regions.  fMRI responses consist in voxels whose three-dimensional
%spatial arrangement respects the anatomy of the brain and correspondingly the discriminative voxels are
%thus expected to have a specific localized spatial organization, which is important
%for the subsequent identification task performed by neuroscientists.  %one is interested in localizing areas in functional magnetic resonance
%imaging (fMRI).
 %


    %matical ignore this spatial configuration, with a
%potential loss in interpretability and performance.




%Recent years have seen a strong growth of interest in multivariate approaches for analysing brain
%activity patterns. The primary goal of these approaches is to reveal the information represented in
%neuronal population codes.  %-information
%analysis

%\comm{Several methods for carrying out high-dimensional inference that are currently in popular use
%rely on marginal measures of the significance of individual variables. An important example
%is the Stability Selection methodology for variable selection proposed in Meinshausen and
%B��uhlmann (2010) and further developed in Shah and Samworth (2013). Such algorithms are
%vulnerable in situations with two or more highly correlated signal variables. In the case of
%Stability Selection, this is because these variables can ��split the vote��, with the result that
%neither/none of the variables is chosen. Similar issues arise with other variable selection
%methods based on p-values for the additional effect of each individual variable. Switching the
%focus to the significance of clusters (determined for instance using the methodology outlined in
%the paper), rather than individual variables, offers the potential t}


%demonstrate the effectiveness of our method

%We first solve the problem of revealing a region��s ��involvement�� in some cognitive function.
%\comm{As a future work, we can consider the sparsity induced norms. }



%\comm{This method can be extended to selection other kind of features beyond voxels, such as network  \cite{Bellec06ClusteringNetwork}.  In more details.  }%



%\section{Proof of the First Zonklar Equation}
%Appendix one text goes here.

% you can choose not to have a title for an appendix
% if you want by leaving the argument blank
%\section{}
%Appendix two text goes here.


% use section* for acknowledgement
\section*{Acknowledgment}
Thanks to Prof. Alexandre Gramfort of Telecom ParisTech  for kindly providing us with  TV-L1 code, which is  under integration in the Nilearn package.

%This work was supported by the Natural Science Foundation of China, Grant
%Nos. 61035006, 61125304, 11201054, 61273361, 91330201, and by  the National Basic Research Program (973 Program), Grant No. 2015CB856000, and by the Specialized Research Fund for the
%Doctoral Program of Higher Education of China 20120185110028, and by the Fundamental Research Funds for the Central Universities ZYGX2012J118, ZYGX2013Z005 and by the Key Technology R\&D Program of Sichuan Province 2012SZ0159.





%% The Appendices part is started with the command \appendix;
%% appendix sections are then done as normal sections
%% \appendix

%% \section{}
%% \label{}

%% References
%%
%% Following citation commands can be used in the body text:
%%
%%  \citet{key}  ==>>  Jones et al. (1990)
%%  \citep{key}  ==>>  (Jones et al., 1990)
%%
%% Multiple citations as normal:
%% \citep{key1,key2}         ==>> (Jones et al., 1990; Smith, 1989)
%%                            or  (Jones et al., 1990, 1991)
%%                            or  (Jones et al., 1990a,b)
%% \cite{key} is the equivalent of \citet{key} in author-year mode
%%
%% Full author lists may be forced with \citet* or \citep*, e.g.
%%   \citep*{key}            ==>> (Jones, Baker, and Williams, 1990)
%%
%% Optional notes as:
%%   \citep[chap. 2]{key}    ==>> (Jones et al., 1990, chap. 2)
%%   \citep[e.g.,][]{key}    ==>> (e.g., Jones et al., 1990)
%%   \citep[see][pg. 34]{key}==>> (see Jones et al., 1990, pg. 34)
%%  (Note: in standard LaTeX, only one note is allowed, after the ref.
%%   Here, one note is like the standard, two make pre- and post-notes.)
%%
%%   \citealt{key}          ==>> Jones et al. 1990
%%   \citealt*{key}         ==>> Jones, Baker, and Williams 1990
%%   \citealp{key}          ==>> Jones et al., 1990
%%   \citealp*{key}         ==>> Jones, Baker, and Williams, 1990
%%
%% Additional citation possibilities
%%   \citeauthor{key}       ==>> Jones et al.
%%   \citeauthor*{key}      ==>> Jones, Baker, and Williams
%%   \citeyear{key}         ==>> 1990
%%   \citeyearpar{key}      ==>> (1990)
%%   \citetext{priv. comm.} ==>> (priv. comm.)
%%   \citenum{key}          ==>> 11 [non-superscripted]
%% Note: full author lists depends on whether the bib style supports them;
%%       if not, the abbreviated list is printed even when full requested.
%%
%% For names like della Robbia at the start of a sentence, use
%%   \Citet{dRob98}         ==>> Della Robbia (1998)
%%   \Citep{dRob98}         ==>> (Della Robbia, 1998)
%%   \Citeauthor{dRob98}    ==>> Della Robbia

\begin{figure*}
    \centering
    \includegraphics[width=0.8\textwidth,height=0.85\textwidth]{T1_compare_9sites_099.eps}
    \vspace{-0.1cm}\caption{Selected overlapped voxels of our method shared by at least 9 centers for the sMRI data, where the probability value is 0.99. We compare the results by different methods. We can see our method is the best which can find out significantly interpretable discriminative structures. Randomized $\ell_1$ method and $\ell_2$ logistic regression are the second best and all the rest methods fail to find out consistent selected voxels among different centers.
   }\label{Fig:Compare9_T1}
\end{figure*}

\begin{figure*}
    \centering
    \includegraphics[width=0.8\textwidth,height=0.85\textwidth]{T1_compare_11sites_099.eps}
    \vspace{-0.1cm}\caption{Selected overlapped voxels of our method shared by at least 11 centers for the sMRI data, where the probability value is 0.99. We compare the results by different methods. We can see our method is the best in terms of revealing   interpretable discriminative structures.
   }\label{Fig:Compare11_T1}
\end{figure*}

\begin{figure*}
    \centering
    \includegraphics[width=0.8\textwidth,height=0.85\textwidth]{T1_compare_14sites_099.eps}
    \vspace{-0.1cm}\caption{Selected overlapped voxels of our method shared by 14 centers for the sMRI data, where the probability value is 0.99. We compare the results by different methods. We can see our method is the  best in terms of revealing   interpretable discriminative structures.
   }\label{Fig:Compare14_T1}
\end{figure*}




\begin{figure*}
    \centering
    \includegraphics[width=0.8\textwidth,height=0.85\textwidth]{zalff_compare_9sites_0975.eps}
    \vspace{-0.1cm}\caption{Selected overlapped voxels of our method shared by at least 9 centers for the fMRI data, where the probability value is 0.975. We compare the results by different methods. We can see our method is the best which can find out significantly interpretable discriminative structures. Randomized $\ell_1$ method and $\ell_2$ logistic regression are the second best and all the rest methods fail to find out consistent selected voxels among different centers.
   }\label{Fig:Compare9_fMRI}
\end{figure*}

\begin{figure*}
    \centering
    \includegraphics[width=0.8\textwidth,height=0.85\textwidth]{zalff_compare_11sites_0975.eps}
    \vspace{-0.1cm}\caption{Selected overlapped voxels of our method shared by at least 11 centers for the fMRI data, where the probability value is 0.975. We compare the results by different methods. We can see our method is the best in terms of revealing   interpretable discriminative structures.
   }\label{Fig:Compare11_fMRI}
\end{figure*}

\begin{figure*}
    \centering
    \includegraphics[width=0.8\textwidth,height=0.85\textwidth]{zalff_compare_14sites_0975.eps}
    \vspace{-0.1cm}\caption{Selected overlapped voxels of our method shared by 14 centers for the fMRI data, where the probability value is 0.975. We compare the results by different methods. We can see our method is the best in terms of revealing   interpretable discriminative structures.
   }\label{Fig:Compare14_fMRI}
\end{figure*}

%% References with bibTeX database:

%\bibliographystyle{model2-names}
\bibliographystyle{IEEEtran}
\bibliography{GroupSparsity,Pattern_recognition,brainScience,StabilitySelection,multi-center}

%% Authors are advised to submit their bibtex database files. They are
%% requested to list a bibtex style file in the manuscript if they do
%% not want to use model2-names.bst.

%% References without bibTeX database:

% \begin{thebibliography}{00}

%% \bibitem must have one of the following forms:
%%   \bibitem[Jones et al.(1990)]{key}...
%%   \bibitem[Jones et al.(1990)Jones, Baker, and Williams]{key}...
%%   \bibitem[Jones et al., 1990]{key}...
%%   \bibitem[\protect\citeauthoryear{Jones, Baker, and Williams}{Jones
%%       et al.}{1990}]{key}...
%%   \bibitem[\protect\citeauthoryear{Jones et al.}{1990}]{key}...
%%   \bibitem[\protect\astroncite{Jones et al.}{1990}]{key}...
%%   \bibitem[\protect\citename{Jones et al., }1990]{key}...
%%   \harvarditem[Jones et al.]{Jones, Baker, and Williams}{1990}{key}...
%%

% \bibitem[ ()]{}

% \end{thebibliography}



\end{document}








\begin{IEEEbiography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]{yilun}}]{Yilun Wang}
earned a doctorate in computational and applied mathematics from Rice University in 2009. From 2009 to 2012, he was a postdoctoral researcher in Cornell University. From 2012 until now, he is an associate professor of school of mathematical sciences of University of Electronic Science and Technology of China (UESTC). Since December 2014,  he is also a research associate of center for applied mathematics of Cornell University.     He  is working on developing %numerical optimization algorithms.
fast algorithms to
find the ``optimal" sparse representation, with applications %in image processing and vision computing, and machine learning. Specifically, I have been working on sparsity applied
in visual computing, machine learning  and data mining to deal with high dimensionality. He has published several papers in high ranked journals such as SIAM Journal on Imaging Sciences, IEEE Transactions on Image Processing, Information Sciences and Journal of agricultural, biological, and environmental statistics. One of his papers has been cited over $780$ times according to Google Scholar.  %, such as sparse, large-scale learning with
\end{IEEEbiography}

% if you will not have a photo at all:
\begin{IEEEbiography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]{zhangsheng}}]{Sheng Zhang}
born in 1989.07, postgraduate student at University of Electronic Science Technology of China (UESTC).
Research areas: Statistical machine learning algorithms; Sparse modeling in neuroimaging; Brain pattern recognition technique and cognitive application; Email: shengzcn@gmail.com
\end{IEEEbiography}

% if you will not have a photo at all:
\begin{IEEEbiography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]{junjie}}]{Junjie Zheng}
born in 1987.08, Ph.d student at University of Electronic Science Technology of China (UESTC) .
Research areas include Brain pattern recognition technique and cognitive application; Dynamic Brain network; Multi- model method and application in Nervous disease and mental disease.
\end{IEEEbiography}

\begin{IEEEbiography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]{Heng}}]{Heng Chen}
born in 1988.01, Ph.d student at University of Electronic Science Technology of China (UESTC) .
Research areas . Deficits of brain network of ASD; Dynamic property of human brain network; Application of MVPA in nervous disease and mental disease. Email: chenheng\_usetc@163.com
\end{IEEEbiography}


\begin{IEEEbiography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]{huafu}}]{Huafu Chen}
%\begin{IEEEbiography}{Michael Shell}
 born in 1967.11, Professor at University of Electronic Science Technology of China(UESTC) . Research areas include Brain pattern recognition technique and cognitive application; Dynamic Brain network; Multi- model method and application in Nervous disease and mental disease. He have published SCI 100 articles including Brain, The American Journal of Psychiatry ,NeuroImage, Human Brian Mapping, PLoS One, IEEE Trans MI/BME etc Journal.
\end{IEEEbiography}

%\end{document}












bare_jrnl.bbl
% Generated by IEEEtran.bst, version: 1.13 (2008/09/30)
\begin{thebibliography}{10}
\providecommand{\url}[1]{#1}
\csname url@samestyle\endcsname
\providecommand{\newblock}{\relax}
\providecommand{\bibinfo}[2]{#2}
\providecommand{\BIBentrySTDinterwordspacing}{\spaceskip=0pt\relax}
\providecommand{\BIBentryALTinterwordstretchfactor}{4}
\providecommand{\BIBentryALTinterwordspacing}{\spaceskip=\fontdimen2\font plus
\BIBentryALTinterwordstretchfactor\fontdimen3\font minus
  \fontdimen4\font\relax}
\providecommand{\BIBforeignlanguage}[2]{{%
\expandafter\ifx\csname l@#1\endcsname\relax
\typeout{** WARNING: IEEEtran.bst: No hyphenation pattern has been}%
\typeout{** loaded for the language `#1'. Using the pattern for}%
\typeout{** the default language instead.}%
\else
\language=\csname l@#1\endcsname
\fi
#2}}
\providecommand{\BIBdecl}{\relax}
\BIBdecl

\bibitem{Haxby01MVA}
J.~V. Haxby, M.~I. Gobbini, M.~L. Furey, A.~Ishai, J.~L. Schouten, and
  P.~Pietrini, ``Distributed and overlapping representations of faces and
  objects in ventral temporal cortex,'' \emph{Science}, vol. 293, no. 5539, pp.
  2425--2430, 2001.

\bibitem{Hilario08biomarker}
\BIBentryALTinterwordspacing
M.~Hilario and A.~Kalousis, ``Approaches to dimensionality reduction in
  proteomic biomarker studies,'' \emph{Briefings in Bioinformatics}, vol.~9,
  no.~2, pp. 102--118, 2008. [Online]. Available:
  \url{http://bib.oxfordjournals.org/content/9/2/102.abstract}
\BIBentrySTDinterwordspacing

\bibitem{Li09SparseVoxelSelection}
Y.~Li, P.~Namburi, Z.~Yu, C.~Guan, J.~Feng, and Z.~Gu, ``Voxel selection in
  fmri data analysis based on sparse representation,'' \emph{Biomedical
  Engineering, IEEE Transactions on}, vol.~56, no.~10, pp. 2439--2451, Oct
  2009.

\bibitem{Guyon03Introduction}
I.~Guyon and A.~Elisseeff, ``An introduction to variable and feature
  selection,'' \emph{Journal of Machine Learning Research}, vol.~3, pp.
  1157--1182, 2003.

\bibitem{Blum1997245}
\BIBentryALTinterwordspacing
A.~L. Blum and P.~Langley, ``Selection of relevant features and examples in
  machine learning,'' \emph{Artificial Intelligence}, vol.~97, no. 1�C2, pp.
  245 -- 271, 1997, relevance. [Online]. Available:
  \url{http://www.sciencedirect.com/science/article/pii/S0004370297000635}
\BIBentrySTDinterwordspacing

\bibitem{Bennett10reliableFMRI}
\BIBentryALTinterwordspacing
C.~M. Bennett and M.~B. Miller, ``How reliable are the results from functional
  magnetic resonance imaging?'' \emph{Annals of the New York Academy of
  Sciences}, vol. 1191, no.~1, pp. 133--155, 2010. [Online]. Available:
  \url{http://dx.doi.org/10.1111/j.1749-6632.2010.05446.x}
\BIBentrySTDinterwordspacing

\bibitem{buhlmann2011statistics}
P.~B{\"u}hlmann and S.~Van De~Geer, \emph{Statistics for high-dimensional data:
  methods, theory and applications}.\hskip 1em plus 0.5em minus 0.4em\relax
  Springer, 2011.

\bibitem{Anderson10critiqueMVPA}
M.~L. Anderson and T.~Oates, ``A critique of multi-voxel pattern analysis,'' in
  \emph{Proceedings of the 32nd Annual Conference of the Cognitive Science
  Society}, 2010.

\bibitem{Yu08denseFeatureGroups}
\BIBentryALTinterwordspacing
L.~Yu, C.~Ding, and S.~Loscalzo, ``Stable feature selection via dense feature
  groups,'' in \emph{Proceedings of the 14th ACM SIGKDD International
  Conference on Knowledge Discovery and Data Mining}, ser. KDD '08.\hskip 1em
  plus 0.5em minus 0.4em\relax New York, NY, USA: ACM, 2008, pp. 803--811.
  [Online]. Available: \url{http://doi.acm.org/10.1145/1401890.1401986}
\BIBentrySTDinterwordspacing

\bibitem{He2010215}
\BIBentryALTinterwordspacing
Z.~He and W.~Yu, ``Stable feature selection for biomarker discovery,''
  \emph{Computational Biology and Chemistry}, vol.~34, no.~4, pp. 215 -- 225,
  2010. [Online]. Available:
  \url{http://www.sciencedirect.com/science/article/pii/S1476927110000502}
\BIBentrySTDinterwordspacing

\bibitem{buhlmann2011controlling}
\BIBentryALTinterwordspacing
P.~B{\"u}hlmann, P.~R{\"u}timann, and M.~Kalisch, ``Controlling false positive
  selections in high-dimensional regression and causal inference,'' \emph{Stat
  Methods Med Res}, vol. online first, Nov. 2011. [Online]. Available:
  \url{http://dx.doi.org/10.1177/0962280211428371}
\BIBentrySTDinterwordspacing

\bibitem{Cover65Pattern}
T.~Cover, ``Geometrical and statistical properties of systems of linear
  inequalities with applications in pattern recognition,'' \emph{IEEE
  transactions on electronic computers}, vol.~14, no.~3, pp. 326--334, 1965.

\bibitem{Raizada2010PFN}
\BIBentryALTinterwordspacing
R.~D.~S. Raizada and N.~Kriegeskorte, ``Pattern-information fmri: New questions
  which it opens up and challenges which face it,'' \emph{Int. J. Imaging Syst.
  Technol.}, vol.~20, no.~1, pp. 31--41, Mar. 2010. [Online]. Available:
  \url{http://dx.doi.org/10.1002/ima.v20:1}
\BIBentrySTDinterwordspacing

\bibitem{Mueller10how}
B.~Mueller, ``How to do a functional multi-center neuroimaging study,'' in
  \emph{Presentation at the Eighteenth Scientific Meeting of the International
  Society for Magnetic Resonance in Medicine}, Stockholm, Sweden, May 2010.

\bibitem{Poldrack14neuroimaging}
R.~A. Poldrack and K.~J. Gorgolewski, ``Making big data open: data sharing in
  neuroimaging,'' \emph{Nat Neurosci}, vol.~17, no.~11, pp. 1510--1517, 2014.

\bibitem{Serres14multimodals}
\BIBentryALTinterwordspacing
S.~Serres, C.~J. Martin, M.~Sarmiento~Soto, C.~Bristow, E.~R. O'Brien, J.~J.
  Connell, A.~A. Khrapitchev, and N.~R. Sibson, ``Structural and functional
  effects of metastases in rat brain determined by multimodal mri,''
  \emph{International Journal of Cancer}, vol. 134, no.~4, pp. 885--896, 2014.
  [Online]. Available: \url{http://dx.doi.org/10.1002/ijc.28406}
\BIBentrySTDinterwordspacing

\bibitem{Wainer14}
J.~Wainer, B.~Robins, F.~Amirabdollahian, and K.~Dautenhahn, ``Using the
  humanoid robot kaspar to autonomously play triadic games and facilitate
  collaborative play among children with autism,'' \emph{Autonomous Mental
  Development, IEEE Transactions on}, vol.~6, no.~3, pp. 183--199, Sept 2014.

\bibitem{Haxby14Neural}
J.~V. Haxby, A.~C. Connolly, and J.~S. Guntupalli, ``Decoding neural
  representational spaces using multivariate pattern analysis,'' \emph{Annual
  Review of Neuroscience}, vol.~37, no.~1, pp. 435--456, 2014.

\bibitem{yu2013}
\BIBentryALTinterwordspacing
B.~Yu, ``Stability,'' \emph{Bernoulli}, vol.~19, no.~4, pp. 1484--1500, 09
  2013. [Online]. Available: \url{http://dx.doi.org/10.3150/13-BEJSP14}
\BIBentrySTDinterwordspacing

\bibitem{Martino13MP}
A.~Di~Martino, C.-G. Yan, Q.~Li, E.~Denio, F.~X. Castellanos, and etc, ``The
  autism brain imaging data exchange: towards a large-scale evaluation of the
  intrinsic brain architecture in autism,'' \emph{Mol Psychiatry}, vol.~19,
  no.~6, pp. 659--667, 2014.

\bibitem{Yamashita08SparseEstimation}
O.~Yamashita, M.~Aki~Sato, T.~Yoshioka, F.~Tong, and Y.~Kamitani, ``Sparse
  estimation automatically selects voxels relevant for the decoding of fmri
  activity patterns,'' \emph{NeuroImage}, vol.~42, no.~4, pp. 1414--1429, 2008.

\bibitem{tibshirani1996regression}
R.~Tibshirani, ``Regression shrinkage and selection via the lasso,''
  \emph{Journal of the Royal Statistical Society. Series B (Methodological)},
  pp. 267--288, 1996.

\bibitem{Rasmussen2012SparsityFMRI}
P.~M. Rasmussen, L.~K. Hansen, K.~H. Madsen, N.~W. Churchill, and S.~C.
  Strother, ``Model sparsity and brain pattern interpretation of classification
  models in neuroimaging,'' \emph{Pattern Recognition}, vol.~45, no.~6, pp.
  2085 -- 2100, 2012.

\bibitem{Tononi94segregation}
G.~Tononi, O.~Sporns, and G.~M. Edelman, ``A measure for brain complexity:
  relating functional segregation and integration in the nervous system,''
  \emph{Proceedings of the National Academy of Sciences}, vol.~91, no.~11, pp.
  5033--5037, 1994.

\bibitem{Xiang:2009}
Z.~Xiang, Y.~Xi, U.~Hasson, and P.~Ramadge, ``{Boosting with Spatial
  Regularization},'' in \emph{Advances in Neural Information Processing Systems
  22}, Y.~Bengio, D.~Schuurmans, J.~Lafferty, C.~K.~I. Williams, and
  A.~Culotta, Eds., 2009, pp. 2107--2115.

\bibitem{Chu02featureselection}
\BIBentryALTinterwordspacing
C.~Chu, A.-L. Hsu, K.-H. Chou, P.~Bandettini, and C.~Lin, ``{Does feature
  selection improve classification accuracy? Impact of sample size and feature
  selection on classification using anatomical magnetic resonance images},''
  \emph{NeuroImage}, vol.~60, no.~1, pp. 59--70, Mar. 2012. [Online].
  Available: \url{http://dx.doi.org/10.1016/j.neuroimage.2011.11.066}
\BIBentrySTDinterwordspacing

\bibitem{LimbrickOldfield20121230}
\BIBentryALTinterwordspacing
E.~H. Limbrick-Oldfield, J.~C. Brooks, R.~J. Wise, F.~Padormo, J.~V. Hajnal,
  C.~F. Beckmann, and M.~A. Ungless, ``Identification and characterisation of
  midbrain nuclei using optimised functional magnetic resonance imaging,''
  \emph{NeuroImage}, vol.~59, no.~2, pp. 1230 -- 1238, 2012. [Online].
  Available:
  \url{http://www.sciencedirect.com/science/article/pii/S1053811911009037}
\BIBentrySTDinterwordspacing

\bibitem{Zou05ElsticNet}
H.~Zou and T.~Hastie, ``Regularization and variable selection via the elastic
  net,'' \emph{Journal of the Royal Statistical Society: Series B (Statistical
  Methodology)}, vol.~67, no.~2, pp. 301--320, 2005.

\bibitem{RyaliCSM12L1L2Network}
S.~Ryali, T.~Chen, K.~Supekar, and V.~Menon, ``Estimation of functional
  connectivity in fmri data using stability selection-based sparse partial
  correlation with elastic net penalty,'' \emph{NeuroImage}, vol.~59, no.~4,
  pp. 3852--3861, 2012.

\bibitem{Ryali10ElsticNetLogistic}
\BIBentryALTinterwordspacing
S.~Ryali, K.~Supekar, D.~A. Abrams, and V.~Menon, ``Sparse logistic regression
  for whole-brain classification of fmri data,'' \emph{NeuroImage}, vol.~51,
  no.~2, pp. 752 -- 764, 2010. [Online]. Available:
  \url{http://www.sciencedirect.com/science/article/pii/S1053811910002089}
\BIBentrySTDinterwordspacing

\bibitem{Gramfort:2013:IPR:2552484.2552506}
\BIBentryALTinterwordspacing
A.~Gramfort, B.~Thirion, and G.~Varoquaux, ``Identifying predictive regions
  from fmri with tv-l1 prior,'' in \emph{Proceedings of the 2013 International
  Workshop on Pattern Recognition in Neuroimaging}, ser. PRNI '13.\hskip 1em
  plus 0.5em minus 0.4em\relax Washington, DC, USA: IEEE Computer Society,
  2013, pp. 17--20. [Online]. Available:
  \url{http://dx.doi.org/10.1109/PRNI.2013.14}
\BIBentrySTDinterwordspacing

\bibitem{Bach12StructuredSparsity}
F.~Bach, R.~Jenatton, J.~Mairal, and G.~Obozinski, ``Structured sparsity
  through convex optimization,'' \emph{Statist. Sci.}, vol.~27, no.~4, p.
  450�C468, 2012.

\bibitem{Baldassarre12StructuredSparsityFMRI}
L.~Baldassarre, J.~Mourao-Miranda, and M.~Pontil, ``Structured sparsity models
  for brain decoding from fmri data,'' in \emph{Pattern Recognition in
  NeuroImaging (PRNI), 2012 International Workshop on}, 2012, pp. 5--8.

\bibitem{Michael11TV}
V.~Michel, A.~Gramfort, G.~Varoquaux, E.~Eger, and B.~Thirion, ``Total
  variation regularization for fmri-based prediction of behavior,''
  \emph{Medical Imaging, IEEE Transactions on}, vol.~30, no.~7, pp. 1328--1340,
  2011.

\bibitem{Ye12groupNonconvex}
S.~Xiang, X.~Shen, and J.~Ye, ``Efficient sparse group feature selection via
  nonconvex optimization,'' \emph{CoRR}, vol. abs/1205.5075, 2012.

\bibitem{Liu10NIPSGroupeTree}
J.~Liu and J.~Ye, ``Moreau-yosida regularization for grouped tree structure
  learning,'' in \emph{Advances in Neural Information Processing Systems 23},
  J.~Lafferty, C.~Williams, J.~Shawe-taylor, R.~Zemel, and A.~Culotta, Eds.,
  2010, pp. 1459--1467.

\bibitem{Yuan13Overlapping}
L.~Yuan, J.~Liu, and J.~Ye, ``Efficient methods for overlapping group lasso,''
  \emph{Pattern Analysis and Machine Intelligence, IEEE Transactions on},
  vol.~PP, no.~99, pp. 1--1, 2013.

\bibitem{Jacob09GLO}
L.~Jacob, G.~Obozinski, and J.-P. Vert, ``Group lasso with overlap and graph
  lasso,'' in \emph{Proceedings of the 26th Annual International Conference on
  Machine Learning}, ser. ICML '09.\hskip 1em plus 0.5em minus 0.4em\relax New
  York, NY, USA: ACM, 2009, pp. 433--440.

\bibitem{Ng11classifier}
B.~Ng and R.~Abugharbieh, ``Generalized group sparse classifiers with
  application in fmri brain decoding,'' in \emph{Proceedings of the 2011 IEEE
  Conference on Computer Vision and Pattern Recognition}, ser. CVPR '11.\hskip
  1em plus 0.5em minus 0.4em\relax Washington, DC, USA: IEEE Computer Society,
  2011, pp. 1065--1071.

\bibitem{Ye09L21}
J.~Liu, S.~Ji, and J.~Ye, ``Multi-task feature learning via efficient
  $\ell_{2,1}$-norm minimization,'' in \emph{Proceedings of the Twenty-Fifth
  Conference on Uncertainty in Artificial Intelligence}, ser. UAI '09.\hskip
  1em plus 0.5em minus 0.4em\relax Arlington, Virginia, United States: AUAI
  Press, 2009, pp. 339--348.

\bibitem{Batmanghelich12Groupsparisty}
N.~Batmanghelich, B.~Taskar, and C.~Davatzikos, ``Generative-discriminative
  basis learning for medical imaging,'' \emph{Medical Imaging, IEEE
  Transactions on}, vol.~31, no.~1, pp. 51--69, Jan 2012.

\bibitem{Michel2012SupervisedClustering}
V.~Michel, A.~Gramfort, G.~Varoquaux, E.~Eger, C.~Keribin, and B.~Thirion, ``A
  supervised clustering approach for fmri-based inference of brain states,''
  \emph{Pattern Recognition}, vol.~45, no.~6, pp. 2041 -- 2049, 2012.

\bibitem{jenatton12multiscale}
R.~Jenatton, A.~Gramfort, V.~Michel, G.~Obozinski, E.~Eger, F.~Bach, and
  B.~Thirion, ``Multiscale mining of fmri data with hierarchical structured
  sparsity,'' \emph{SIAM J. Imaging Sciences}, pp. 835--856, 2012.

\bibitem{Meinshausen10StabilitySelection}
\BIBentryALTinterwordspacing
N.~Meinshausen and P.~B��hlmann, ``Stability selection,'' \emph{Journal of the
  Royal Statistical Society: Series B (Statistical Methodology)}, vol.~72,
  no.~4, pp. 417--473, 2010. [Online]. Available:
  \url{http://dx.doi.org/10.1111/j.1467-9868.2010.00740.x}
\BIBentrySTDinterwordspacing

\bibitem{Shah13StabilitySelection}
R.~D. Shah and R.~J. Samworth, ``Variable selection with error control: another
  look at stability selection,'' \emph{J. R. Statist. Soc. B}, vol.~75, no.
  part 1, pp. 55--80, 2013.

\bibitem{Hastie09Ensemle}
T.~Hastie, R.~Tibshirani, and J.~Friedman,
  ``\BIBforeignlanguage{English}{Ensemble learning},'' in
  \emph{\BIBforeignlanguage{English}{The Elements of Statistical Learning}},
  ser. Springer Series in Statistics.\hskip 1em plus 0.5em minus 0.4em\relax
  Springer New York, 2009, pp. 605--624.

\bibitem{Mitchell04LearningDecodeCognitive}
T.~M. Mitchell, R.~Hutchinson, R.~S. Niculescu, F.~Pereira, X.~Wang, M.~Just,
  and S.~Newman, ``Learning to decode cognitive states from brain images,''
  \emph{Mach. Learn.}, vol.~57, no. 1-2, pp. 145--175, Oct. 2004.

\bibitem{Li12SparseLocalization}
Y.~Li, J.~Long, L.~He, H.~Lu, Z.~Gu, and et~al, ``A sparse representation-based
  algorithm for pattern localization in brain imaging data analysis,''
  \emph{PLoS ONE}, vol.~7, no.~12, 2012.

\bibitem{Poldrack06Criti}
R.~A. Poldrack, ``Can cognitive processes be inferred from neuroimaging data?''
  \emph{Trends in cognitive sciences}, vol.~10, no.~2, p.~59, 63 2006.

\bibitem{Ye12SparseLearningSS}
J.~Ye, M.~Farnum, E.~Yang, R.~Verbeeck, V.~Lobanov, N.~Raghavan, G.~Novak,
  A.~DiBernardo, V.~Narayan, and for~the Alzheimer's Disease
  Neuroimaging~Initiative, ``Sparse learning and stability selection for
  predicting mci to ad conversion using baseline adni data,'' \emph{BMC
  Neurology}, vol.~12, no.~1, p.~46, 2012.

\bibitem{Cao2014Sparse}
H.~Cao, J.~Duan, D.~Lin, Y.~Y. Shugart, V.~Calhoun, and Y.-P. Wang, ``Sparse
  representation-based biomarker selection for schizophrenia with integrated
  analysis of fmri and \{SNPs\},'' \emph{NeuroImage}, no.~0, pp.~--, 2014.

\bibitem{Ryali20123852}
\BIBentryALTinterwordspacing
S.~Ryali, T.~Chen, K.~Supekar, and V.~Menon, ``Estimation of functional
  connectivity in fmri data using stability selection-based sparse partial
  correlation with elastic net penalty,'' \emph{NeuroImage}, vol.~59, no.~4,
  pp. 3852 -- 3861, 2012. [Online]. Available:
  \url{http://www.sciencedirect.com/science/article/pii/S105381191101336X}
\BIBentrySTDinterwordspacing

\bibitem{Wang14RSS}
Y.~Wang, J.~Zheng, S.~Zhang, X.~Duan, and H.~Chen, ``Randomized structural
  sparsity via constrained block subsampling for improved sensitivity of
  discriminative voxel identification,'' \emph{arXiv}, 2014.

\bibitem{Li13ClusterGuided}
Z.~Li, J.~Liu, Y.~Yang, X.~Zhou, and H.~Lu, ``Clustering-guided sparse
  structural learning for unsupervised feature selection,'' \emph{IEEE
  Transactions on Knowledge and Data Engineering}, vol.~99, no. PrePrints,
  p.~1, 2013.

\bibitem{Mairal13SFS}
\BIBentryALTinterwordspacing
J.~Mairal and B.~Yu, ``Supervised feature selection in graphs with path coding
  penalties and network flows,'' \emph{J. Mach. Learn. Res.}, vol.~14, no.~1,
  pp. 2449--2485, Jan. 2013. [Online]. Available:
  \url{http://dl.acm.org/citation.cfm?id=2567709.2567740}
\BIBentrySTDinterwordspacing

\bibitem{Bach12OSP}
F.~Bach, R.~Jenatton, J.~Mairal, and G.~Obozinski, ``Optimization with
  sparsity-inducing penalties,'' \emph{Found. Trends Mach. Learn.}, vol.~4,
  no.~1, pp. 1--106, 2012.

\bibitem{Beinrucker12ExtensionSS}
A.~Beinrucker, u.~Dogan, and G.~Blanchard, ``A simple extension of stability
  feature selection,'' in \emph{Pattern Recognition}, ser. Lecture Notes in
  Computer Science, A.~Pinz, T.~Pock, H.~Bischof, and F.~Leberl, Eds.\hskip 1em
  plus 0.5em minus 0.4em\relax Springer Berlin Heidelberg, 2012, vol. 7476, pp.
  256--265.

\bibitem{Lahiri99Blockbootstrap}
S.~N. Lahiri, ``Theoretical comparisons of block bootstrap methods,''
  \emph{Ann. Statist.}, vol.~27, no.~1, pp. 386--404, 1999.

\bibitem{TzourioMazoyer2002AnatomicalParcellation}
\BIBentryALTinterwordspacing
N.~Tzourio-Mazoyer, B.~Landeau, D.~Papathanassiou, F.~Crivello, O.~Etard,
  N.~Delcroix, B.~Mazoyer, and M.~Joliot, ``Automated anatomical labeling of
  activations in \{SPM\} using a macroscopic anatomical parcellation of the
  \{MNI\} \{MRI\} single-subject brain,'' \emph{NeuroImage}, vol.~15, no.~1,
  pp. 273 -- 289, 2002. [Online]. Available:
  \url{http://www.sciencedirect.com/science/article/pii/S1053811901909784}
\BIBentrySTDinterwordspacing

\bibitem{Sarndal03stratified}
C.-E. e.~a. S\"{a}rndal, \emph{Model Assisted Survey Sampling}.\hskip 1em plus
  0.5em minus 0.4em\relax Springer, 2003, ch. Stratified Sampling, pp.
  100--109.

\bibitem{Devries86}
P.~de~Vries, ``\BIBforeignlanguage{English}{Stratified random sampling},'' in
  \emph{\BIBforeignlanguage{English}{Sampling Theory for Forest
  Inventory}}.\hskip 1em plus 0.5em minus 0.4em\relax Springer Berlin
  Heidelberg, 1986, pp. 31--55.

\bibitem{Friedman2006varability}
\BIBentryALTinterwordspacing
L.~Friedman, G.~H. Glover, D.~Krenz, and V.~Magnotta, ``Reducing inter-scanner
  variability of activation in a multicenter fmri study: Role of smoothness
  equalization,'' \emph{NeuroImage}, vol.~32, no.~4, pp. 1656 -- 1668, 2006.
  [Online]. Available:
  \url{http://www.sciencedirect.com/science/article/pii/S1053811906004435}
\BIBentrySTDinterwordspacing

\bibitem{Friedman06JMRI}
\BIBentryALTinterwordspacing
L.~Friedman and G.~H. Glover, ``Report on a multicenter fmri quality assurance
  protocol,'' \emph{Journal of Magnetic Resonance Imaging}, vol.~23, no.~6, pp.
  827--839, 2006. [Online]. Available:
  \url{http://dx.doi.org/10.1002/jmri.20583}
\BIBentrySTDinterwordspacing

\bibitem{Geman:1992:NNB}
\BIBentryALTinterwordspacing
S.~Geman, E.~Bienenstock, and R.~Doursat, ``Neural networks and the
  bias/variance dilemma,'' \emph{Neural Comput.}, vol.~4, no.~1, pp. 1--58,
  Jan. 1992. [Online]. Available:
  \url{http://dx.doi.org/10.1162/neco.1992.4.1.1}
\BIBentrySTDinterwordspacing

\bibitem{James13StasiticalLearning}
G.~James, D.~Witten, T.~Hastie, and R.~Tibshirani, \emph{An Introduction to
  Statistical Learning with Applications in R}, ser. Springer Texts in
  Statistics.\hskip 1em plus 0.5em minus 0.4em\relax Springer, 2013, vol. 103.

\bibitem{Varoquaux12Clustering}
B.~T. Ga\"{e}l~Varoquaux, Alexandre~Gramfort, ``Small-sample brain mapping:
  sparse recovery on spatially correlated designs with randomization and
  clustering,'' in \emph{ICML}, 2012.

\bibitem{Park01042007}
M.~Y. Park, T.~Hastie, and R.~Tibshirani, ``Averaged gene expressions for
  regression,'' \emph{Biostatistics}, vol.~8, no.~2, pp. 212--227, 2007.

\bibitem{REF08a}
R.-E. Fan, K.-W. Chang, C.-J. Hsieh, X.-R. Wang, and C.-J. Lin, ``{LIBLINEAR}:
  A library for large linear classification,'' \emph{Journal of Machine
  Learning Research}, vol.~9, pp. 1871--1874, 2008.

\bibitem{Liu09:SLEP:manual}
\BIBentryALTinterwordspacing
J.~Liu, S.~Ji, and J.~Ye, \emph{SLEP: Sparse Learning with Efficient
  Projections}, Arizona State University, 2009. [Online]. Available:
  \url{http://www.public.asu.edu/~jye02/Software/SLEP}
\BIBentrySTDinterwordspacing

\bibitem{YuFeng2007ALFF}
\BIBentryALTinterwordspacing
Z.~Yu-Feng, H.~Yong, Z.~Chao-Zhe, C.~Qing-Jiu, S.~Man-Qiu, L.~Meng, T.~Li-Xia,
  J.~Tian-Zi, and W.~Yu-Feng, ``Altered baseline brain activity in children
  with \{ADHD\} revealed by resting-state functional \{MRI\},'' \emph{Brain and
  Development}, vol.~29, no.~2, pp. 83 -- 91, 2007. [Online]. Available:
  \url{http://www.sciencedirect.com/science/article/pii/S0387760406001549}
\BIBentrySTDinterwordspacing

\bibitem{Ashburner20008VBM}
\BIBentryALTinterwordspacing
J.~Ashburner and K.~J. Friston, ``Voxel-based morphometry��the methods,''
  \emph{NeuroImage}, vol.~11, no.~6, pp. 805 -- 821, 2000. [Online]. Available:
  \url{http://www.sciencedirect.com/science/article/pii/S1053811900905822}
\BIBentrySTDinterwordspacing

\bibitem{Rondina14}
J.~Rondina, T.~Hahn, L.~de~Oliveira, A.~Marquand, T.~Dresler, T.~Leitner,
  A.~Fallgatter, J.~Shawe-Taylor, and J.~Mourao-Miranda, ``Scors - a method
  based on stability for feature selection and mapping in neuroimaging,''
  \emph{Medical Imaging, IEEE Transactions on}, vol.~33, no.~1, pp. 85--98, Jan
  2014.

\bibitem{Friedman08retest}
\BIBentryALTinterwordspacing
L.~Friedman, H.~Stern, G.~G. Brown, D.~H. Mathalon, J.~Turner, G.~H. Glover,
  R.~L. Gollub, J.~Lauriello, K.~O. Lim, T.~Cannon, D.~N. Greve, H.~J.
  Bockholt, A.~Belger, B.~Mueller, M.~J. Doty, J.~He, W.~Wells, P.~Smyth,
  S.~Pieper, S.~Kim, M.~Kubicki, M.~Vangel, and S.~G. Potkin, ``Test�cretest
  and between-site reliability in a multicenter fmri study,'' \emph{Human Brain
  Mapping}, vol.~29, no.~8, pp. 958--972, 2008. [Online]. Available:
  \url{http://dx.doi.org/10.1002/hbm.20440}
\BIBentrySTDinterwordspacing

\bibitem{wang2010sparse}
Y.~Wang and W.~Yin, ``Sparse signal reconstruction via iterative support
  detection,'' \emph{SIAM Journal on Imaging Sciences}, vol.~3, no.~3, pp.
  462--491, 2010.

\bibitem{Sridharan26082008}
\BIBentryALTinterwordspacing
D.~Sridharan, D.~J. Levitin, and V.~Menon, ``A critical role for the right
  fronto-insular cortex in switching between central-executive and default-mode
  networks,'' \emph{Proceedings of the National Academy of Sciences}, vol. 105,
  no.~34, pp. 12\,569--12\,574, 2008. [Online]. Available:
  \url{http://www.pnas.org/content/105/34/12569.abstract}
\BIBentrySTDinterwordspacing

\bibitem{Bonnelle20032012}
\BIBentryALTinterwordspacing
V.~Bonnelle, T.~E. Ham, R.~Leech, K.~M. Kinnunen, M.~A. Mehta, R.~J. Greenwood,
  and D.~J. Sharp, ``Salience network integrity predicts default mode network
  function after traumatic brain injury,'' \emph{Proceedings of the National
  Academy of Sciences}, vol. 109, no.~12, pp. 4690--4695, 2012. [Online].
  Available: \url{http://www.pnas.org/content/109/12/4690.abstract}
\BIBentrySTDinterwordspacing

\bibitem{Schumann14072004}
\BIBentryALTinterwordspacing
C.~M. Schumann, J.~Hamstra, B.~L. Goodlin-Jones, L.~J. Lotspeich, H.~Kwon,
  M.~H. Buonocore, C.~R. Lammers, A.~L. Reiss, and D.~G. Amaral, ``The amygdala
  is enlarged in children but not adolescents with autism; the hippocampus is
  enlarged at all ages,'' \emph{The Journal of Neuroscience}, vol.~24, no.~28,
  pp. 6392--6401, 2004. [Online]. Available:
  \url{http://www.jneurosci.org/content/24/28/6392.abstract}
\BIBentrySTDinterwordspacing

\bibitem{Gusnard27032001}
\BIBentryALTinterwordspacing
D.~A. Gusnard, E.~Akbudak, G.~L. Shulman, and M.~E. Raichle, ``Medial
  prefrontal cortex and self-referential mental activity: Relation to a default
  mode of brain function,'' \emph{Proceedings of the National Academy of
  Sciences}, vol.~98, no.~7, pp. 4259--4264, 2001. [Online]. Available:
  \url{http://www.pnas.org/content/98/7/4259.abstract}
\BIBentrySTDinterwordspacing

\bibitem{Bush2000215}
\BIBentryALTinterwordspacing
G.~Bush, P.~Luu, and M.~I. Posner, ``Cognitive and emotional influences in
  anterior cingulate cortex,'' \emph{Trends in Cognitive Sciences}, vol.~4,
  no.~6, pp. 215 -- 222, 2000. [Online]. Available:
  \url{http://www.sciencedirect.com/science/article/pii/S1364661300014832}
\BIBentrySTDinterwordspacing

\bibitem{Castelli02}
F.~Castelli, C.~Frith, F.~Happ{\'e}, and U.~Frith, ``Autism, asperger syndrome
  and brain mechanisms for the attribution of mental states to animated
  shapes,'' \emph{Brain}, vol. 125, no.~8, pp. 1839--1849, 2002.

\bibitem{Critchley00}
H.~D. Critchley, E.~M. Daly, E.~T. Bullmore, S.~C.~R. Williams,
  T.~Van~Amelsvoort, D.~M. Robertson, A.~Rowe, M.~Phillips, G.~McAlonan,
  P.~Howlin, and D.~G.~M. Murphy, ``The functional neuroanatomy of social
  behaviour,'' \emph{Brain}, vol. 123, no.~11, pp. 2203--2212, 2000.

\end{thebibliography}
