Dynamic Data-Driven Estimation of Non-Parametric Choice Models
We study non-parametric estimation of consumer choice models. Non-parametric approaches were introduced to alleviate unreasonable assumptions and issues of suboptimal model fit/selection present in traditional parametric approaches, and are prevalent in several application areas. We present two convex optimization-based frameworks, a primal approach and a dual approach, to efficiently learn a non-parametric choice model from data that is close to the best-fitting one. As opposed to the existing literature, both approaches enjoy provable convergence guarantees and extend naturally to the dynamic observation setting. Our computational study on the dynamic setting reveals the true impact of how much data are needed and at what rate to achieve the best trade-off in terms of estimation accuracy and model simplicity. In the static setup, we also compare our non-parametric approach with existing parametric approaches.
comparefig1.tex
\newcommand{\barwid}{.1cm}
\newcommand{\wid}{0.3}
\newcommand{\hgt}{0.35}
\newcommand{\scalefactor}{1}

\pgfplotsset{scaled y ticks=false}
\pgfplotsset{
	/pgfplots/ybar legend/.style={
    /pgfplots/legend image code/.code={%
       \draw[##1,/tikz/.cd,yshift=-0.25em]
        (0cm,0cm) rectangle (3pt,0.8em);},
   },
}
\begin{figure}[t!b]
\centering
\begin{tikzpicture}[scale=\scalefactor]
\begin{axis}[
name=K1,
ybar,
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
symbolic x coords={10,20,50},
xtick=data,
yticklabel style={/pgf/number format/fixed},
ymin=0,ymax=0.07,
xlabel={$m$},
ylabel={Test set MAE},
title={$K=1$},
enlarge x limits=0.25,
]
\addplot table [x=train_m, y=1_nonparametric, col sep=comma] {MAE_test_MMNL_n10_L5_K1.csv};
\addplot table [x=train_m, y=2_nonparametric, col sep=comma] {MAE_test_MMNL_n10_L5_K1.csv};
\addplot table [x=train_m, y=inf_nonparametric, col sep=comma] {MAE_test_MMNL_n10_L5_K1.csv};
\end{axis}

\begin{axis}[
name=K5,
at={(0.25\linewidth, 0)},
ybar,
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
symbolic x coords={10,20,50},
xtick=data,
yticklabel style={/pgf/number format/fixed},
ymajorticks=false,
ymin=0,ymax=0.05,
xlabel={$m$},
legend style={at={(0.5,1.03)},
			  anchor=south,
			  legend columns=3},
title={$K=5$},
enlarge x limits=0.25,
]
\addplot table [x=train_m, y=1_nonparametric, col sep=comma] {MAE_test_MMNL_n10_L5_K5.csv};
\addplot table [x=train_m, y=2_nonparametric, col sep=comma] {MAE_test_MMNL_n10_L5_K5.csv};
\addplot table [x=train_m, y=inf_nonparametric, col sep=comma] {MAE_test_MMNL_n10_L5_K5.csv};
\end{axis}

\begin{axis}[
name=K10,
at={(0.5\linewidth, 0)},
ybar,
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
symbolic x coords={10,20,50},
xtick=data,
yticklabel style={/pgf/number format/fixed},
ymajorticks=false,
ymin=0,ymax=0.05,
xlabel={$m$},
legend style={legend pos=outer north east,
			  cells={anchor=west}},
title={$K=10$},
enlarge x limits=0.25,
]
\addplot table [x=train_m, y=1_nonparametric, col sep=comma] {MAE_test_MMNL_n10_L5_K10.csv};
\addplot table [x=train_m, y=2_nonparametric, col sep=comma] {MAE_test_MMNL_n10_L5_K10.csv};
\addplot table [x=train_m, y=inf_nonparametric, col sep=comma] {MAE_test_MMNL_n10_L5_K10.csv};
\legend{$\ell_1$, $\ell_2$, $\ell_\infty$}
\end{axis}
\end{tikzpicture}
\caption{Test set MAE for $L=5$.}
\label{fig:MAE-m_vs_p}
\end{figure}

comparefig2.tex
\renewcommand{\barwid}{0.15cm}
\renewcommand{\wid}{0.21}
\renewcommand{\hgt}{0.21}
\renewcommand{\scalefactor}{1}
\newcommand{\MAEexpand}{0.001}
\newcommand{\numrankexpand}{100}
\newcommand{\barshift}{3pt}

\pgfplotsset{every tick label/.append style={font=\footnotesize}}
\begin{figure}[h!]
\captionsetup[subfigure]{aboveskip=100pt,belowskip=-1pt}
\centering
\subfigure[$L=5,K=1$]{
\begin{tikzpicture}[scale=\scalefactor]
\pgfplotsset{set layers}
\begin{axis}[
name=L5K1,
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
xticklabels={$\ell_1$,$\ell_2$,$\ell_\infty$},
enlarge x limits=0.25,
y tick label style={
    /pgf/number format/.cd,
        fixed,
        fixed zerofill,
        precision=3,
    /tikz/.cd
},
axis y line*=left,
enlarge y limits={abs=\MAEexpand},
ylabel={\footnotesize Test MAE},
]
\addplot[ybar, ybar legend, /pgf/bar shift=-\barshift, fill=blue!40] table [x=plabel, y=MAE_test, col sep=comma] {nonpar_MAE_test_m20_MMNL_n10_L5_K1.csv};
\label{testMAE}
\end{axis}
\begin{axis}[
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
enlarge x limits=0.25,
xmajorticks=false,
yticklabel style={/pgf/number format/fixed},
axis y line*=right,
ylabel style={align=center},
enlarge y limits={abs=\numrankexpand},
]
\addplot[ybar, ybar legend, /pgf/bar shift=\barshift, fill=red!40] table [x=plabel, y=num_rankings, col sep=comma] {nonpar_MAE_test_m20_MMNL_n10_L5_K1.csv};
\label{numrankings}
\end{axis}
\end{tikzpicture}
\label{fig:compare_L5K1}
}%
\hfill
\subfigure[$L=5,K=5$]{
\begin{tikzpicture}[scale=\scalefactor]
\begin{axis}[
name=L5K5,
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
xticklabels={$\ell_1$,$\ell_2$,$\ell_\infty$},
enlarge x limits=0.25,
y tick label style={
    /pgf/number format/.cd,
        fixed,
        fixed zerofill,
        precision=3,
    /tikz/.cd
},
axis y line*=left,
enlarge y limits={abs=\MAEexpand},
]
\addplot[ybar, ybar legend, /pgf/bar shift=-\barshift, fill=blue!40] table [x=plabel, y=MAE_test, col sep=comma] {nonpar_MAE_test_m20_MMNL_n10_L5_K5.csv};
\end{axis}
\begin{axis}[
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
enlarge x limits=0.25,
xmajorticks=false,
yticklabel style={/pgf/number format/fixed},
axis y line*=right,
ylabel style={align=center},
enlarge y limits={abs=\numrankexpand},
]
\addplot[ybar, ybar legend, /pgf/bar shift=\barshift, fill=red!40] table [x=plabel, y=num_rankings, col sep=comma] {nonpar_MAE_test_m20_MMNL_n10_L5_K5.csv};
\end{axis}
\end{tikzpicture}
\label{fig:compare_L5K5}
}%
\hfill
\subfigure[$L=5,K=10$]{
\begin{tikzpicture}[scale=\scalefactor]
\begin{axis}[
name=L5K10,
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
xticklabels={$\ell_1$,$\ell_2$,$\ell_\infty$},
enlarge x limits=0.25,
y tick label style={
    /pgf/number format/.cd,
        fixed,
        fixed zerofill,
        precision=3,
    /tikz/.cd
},
axis y line*=left,
enlarge y limits={abs=\MAEexpand},
]
\addplot[ybar, /pgf/bar shift=-\barshift, fill=blue!40] table [x=plabel, y=MAE_test, col sep=comma] {nonpar_MAE_test_m20_MMNL_n10_L5_K10.csv};
\end{axis}
\begin{axis}[
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
enlarge x limits=0.25,
xmajorticks=false,
yticklabel style={/pgf/number format/fixed},
axis y line*=right,
ylabel style={align=center},
enlarge y limits={abs=\numrankexpand},
ylabel={\footnotesize Avg \#rankings},
]
\addplot[ybar, /pgf/bar shift=\barshift, fill=red!40] table [x=plabel, y=num_rankings, col sep=comma] {nonpar_MAE_test_m20_MMNL_n10_L5_K10.csv};
\end{axis}
\end{tikzpicture}
\label{fig:compare_L5K10}
}%
\hfill
\subfigure[$L=10,K=1$]{
\begin{tikzpicture}[scale=\scalefactor]
\begin{axis}[
name=L5K5,
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
xticklabels={$\ell_1$,$\ell_2$,$\ell_\infty$},
enlarge x limits=0.25,
y tick label style={
    /pgf/number format/.cd,
        fixed,
        fixed zerofill,
        precision=3,
    /tikz/.cd
},
axis y line*=left,
enlarge y limits={abs=\MAEexpand},
ylabel={\footnotesize Test MAE},
]
\addplot[ybar, /pgf/bar shift=-\barshift, fill=blue!40] table [x=plabel, y=MAE_test, col sep=comma] {nonpar_MAE_test_m20_MMNL_n10_L10_K1.csv};
\end{axis}
\begin{axis}[
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
enlarge x limits=0.25,
xmajorticks=false,
yticklabel style={/pgf/number format/fixed},
axis y line*=right,
ylabel style={align=center},
enlarge y limits={abs=\numrankexpand},
]
\addplot[ybar, /pgf/bar shift=\barshift, fill=red!40] table [x=plabel, y=num_rankings, col sep=comma] {nonpar_MAE_test_m20_MMNL_n10_L10_K1.csv};
\end{axis}
\end{tikzpicture}
\label{fig:compare_L10K1}
}%
\hfill
\subfigure[$L=10,K=5$]{
\begin{tikzpicture}[scale=\scalefactor]
\begin{axis}[
name=L10K5,
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
xticklabels={$\ell_1$,$\ell_2$,$\ell_\infty$},
enlarge x limits=0.25,
y tick label style={
    /pgf/number format/.cd,
        fixed,
        fixed zerofill,
        precision=3,
    /tikz/.cd
},
axis y line*=left,
enlarge y limits={abs=\MAEexpand},
]
\addplot[ybar, /pgf/bar shift=-\barshift, fill=blue!40] table [x=plabel, y=MAE_test, col sep=comma] {nonpar_MAE_test_m20_MMNL_n10_L10_K5.csv};
\end{axis}
\begin{axis}[
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
enlarge x limits=0.25,
xmajorticks=false,
yticklabel style={/pgf/number format/fixed},
axis y line*=right,
ylabel style={align=center},
enlarge y limits={abs=\numrankexpand},
]
\addplot[ybar, /pgf/bar shift=\barshift, fill=red!40] table [x=plabel, y=num_rankings, col sep=comma] {nonpar_MAE_test_m20_MMNL_n10_L10_K5.csv};
\end{axis}
\end{tikzpicture}
\label{fig:compare_L10K5}
}%
\hfill
\subfigure[$L=10,K=10$]{
\begin{tikzpicture}[scale=\scalefactor]
\begin{axis}[
name=L5K5,
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
xticklabels={$\ell_1$,$\ell_2$,$\ell_\infty$},
enlarge x limits=0.25,
y tick label style={
    /pgf/number format/.cd,
        fixed,
        fixed zerofill,
        precision=3,
    /tikz/.cd
},
axis y line*=left,
enlarge y limits={abs=\MAEexpand},
]
\addplot[ybar, /pgf/bar shift=-\barshift, fill=blue!40] table [x=plabel, y=MAE_test, col sep=comma] {nonpar_MAE_test_m20_MMNL_n10_L10_K10.csv};
\end{axis}
\begin{axis}[
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
enlarge x limits=0.25,
xmajorticks=false,
yticklabel style={/pgf/number format/fixed},
axis y line*=right,
ylabel style={align=center},
enlarge y limits={abs=\numrankexpand},
ylabel={\footnotesize Avg \#rankings},
]
\addplot[ybar, /pgf/bar shift=\barshift, fill=red!40] table [x=plabel, y=num_rankings, col sep=comma] {nonpar_MAE_test_m20_MMNL_n10_L10_K10.csv};
\end{axis}
\end{tikzpicture}
\label{fig:compare_L10K10}
}%
\hfill
\subfigure[$L=100,K=1$]{
\begin{tikzpicture}[scale=\scalefactor]
\begin{axis}[
name=L5K5,
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
xticklabels={$\ell_1$,$\ell_2$,$\ell_\infty$},
enlarge x limits=0.25,
y tick label style={
    /pgf/number format/.cd,
        fixed,
        fixed zerofill,
        precision=3,
    /tikz/.cd
},
axis y line*=left,
enlarge y limits={abs=\MAEexpand},
ylabel={\footnotesize Test MAE},
]
\addplot[ybar, /pgf/bar shift=-\barshift, fill=blue!40] table [x=plabel, y=MAE_test, col sep=comma] {nonpar_MAE_test_m20_MMNL_n10_L100_K1.csv};
\end{axis}
\begin{axis}[
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
enlarge x limits=0.25,
xmajorticks=false,
yticklabel style={/pgf/number format/fixed},
axis y line*=right,
ylabel style={align=center},
enlarge y limits={abs=\numrankexpand},
]
\addplot[ybar, /pgf/bar shift=\barshift, fill=red!40] table [x=plabel, y=num_rankings, col sep=comma] {nonpar_MAE_test_m20_MMNL_n10_L100_K1.csv};
\end{axis}
\end{tikzpicture}
\label{fig:compare_L100K1}
}%
\hfill
\subfigure[$L=100,K=5$]{
\begin{tikzpicture}[scale=\scalefactor]
\begin{axis}[
name=L5K5,
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
xticklabels={$\ell_1$,$\ell_2$,$\ell_\infty$},
enlarge x limits=0.25,
y tick label style={
    /pgf/number format/.cd,
        fixed,
        fixed zerofill,
        precision=3,
    /tikz/.cd
},
axis y line*=left,
enlarge y limits={abs=\MAEexpand},
]
\addplot[ybar, /pgf/bar shift=-\barshift, fill=blue!40] table [x=plabel, y=MAE_test, col sep=comma] {nonpar_MAE_test_m20_MMNL_n10_L100_K5.csv};
\end{axis}
\begin{axis}[
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
enlarge x limits=0.25,
xmajorticks=false,
yticklabel style={/pgf/number format/fixed},
axis y line*=right,
ylabel style={align=center},
enlarge y limits={abs=\numrankexpand},
]
\addplot[ybar, /pgf/bar shift=\barshift, fill=red!40] table [x=plabel, y=num_rankings, col sep=comma] {nonpar_MAE_test_m20_MMNL_n10_L100_K5.csv};
\end{axis}
\end{tikzpicture}
\label{fig:compare_L100K5}
}%
\hfill
\subfigure[$L=100,K=10$]{
\begin{tikzpicture}[scale=\scalefactor]
\begin{axis}[
name=L5K5,
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
xticklabels={$\ell_1$,$\ell_2$,$\ell_\infty$},
enlarge x limits=0.25,
y tick label style={
    /pgf/number format/.cd,
        fixed,
        fixed zerofill,
        precision=3,
    /tikz/.cd
},
axis y line*=left,
enlarge y limits={abs=\MAEexpand},
]
\addplot[ybar, /pgf/bar shift=-\barshift, fill=blue!40] table [x=plabel, y=MAE_test, col sep=comma] {nonpar_MAE_test_m20_MMNL_n10_L100_K10.csv};
\end{axis}
\begin{axis}[
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
enlarge x limits=0.25,
xmajorticks=false,
yticklabel style={/pgf/number format/fixed},
axis y line*=right,
ylabel style={align=center},
enlarge y limits={abs=\numrankexpand},
ylabel={\footnotesize Avg \#rankings},
]
\addplot[ybar, /pgf/bar shift=\barshift, fill=red!40] table [x=plabel, y=num_rankings, col sep=comma] {nonpar_MAE_test_m20_MMNL_n10_L100_K10.csv};
\end{axis}
\end{tikzpicture}
\label{fig:compare_L100K10}
}%
\caption{\ref{testMAE} Test set MAE (left axis);\ \ref{numrankings} Avg. \#rankings (right axis).}
\label{fig:MAE_numrankings_compare}
\end{figure}


comparefig3.tex
\renewcommand{\barwid}{0.1cm}
\renewcommand{\wid}{0.45}
\renewcommand{\hgt}{0.35}
\renewcommand{\scalefactor}{1}

\begin{figure}[htb]
\centering
\begin{tikzpicture}[scale=\scalefactor]
\begin{axis}[
name=m10,
ybar,
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
symbolic x coords={1,5,10},
xtick=data,
yticklabel style={/pgf/number format/fixed},
ymin=0,ymax=0.05,
xlabel={$K$},
ylabel={Test set MAE},
title={$m=10$},
enlarge x limits=0.25,
]
\addplot table [x=MMNL_segments_true, y=1_nonparametric, col sep=comma] {MAE_test_m10_MMNL_n10_L5.csv};
\addplot table [x=MMNL_segments_true, y=2_nonparametric, col sep=comma] {MAE_test_m10_MMNL_n10_L5.csv};
\addplot table [x=MMNL_segments_true, y=inf_nonparametric, col sep=comma] {MAE_test_m10_MMNL_n10_L5.csv};
\addplot table [x=MMNL_segments_true, y=1_parametric, col sep=comma] {MAE_test_m10_MMNL_n10_L5.csv};
\addplot table [x=MMNL_segments_true, y=3_parametric, col sep=comma] {MAE_test_m10_MMNL_n10_L5.csv};
\addplot table [x=MMNL_segments_true, y=5_parametric, col sep=comma] {MAE_test_m10_MMNL_n10_L5.csv};
\addplot table [x=MMNL_segments_true, y=10_parametric, col sep=comma] {MAE_test_m10_MMNL_n10_L5.csv};
\end{axis}

\begin{axis}[
name=m20,
at={(0.45\linewidth, 0)},
ybar,
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
symbolic x coords={1,5,10},
xtick=data,
yticklabel style={/pgf/number format/fixed},
ymajorticks=false,
ymin=0,ymax=0.05,
xlabel={$K$},
legend style={at={(0.5,1.03)},
			  anchor=south,
			  legend columns=3},
title={$m=20$},
enlarge x limits=0.25,
]
\addplot table [x=MMNL_segments_true, y=1_nonparametric, col sep=comma] {MAE_test_m20_MMNL_n10_L5.csv};
\addplot table [x=MMNL_segments_true, y=2_nonparametric, col sep=comma] {MAE_test_m20_MMNL_n10_L5.csv};
\addplot table [x=MMNL_segments_true, y=inf_nonparametric, col sep=comma] {MAE_test_m20_MMNL_n10_L5.csv};
\addplot table [x=MMNL_segments_true, y=1_parametric, col sep=comma] {MAE_test_m20_MMNL_n10_L5.csv};
\addplot table [x=MMNL_segments_true, y=3_parametric, col sep=comma] {MAE_test_m20_MMNL_n10_L5.csv};
\addplot table [x=MMNL_segments_true, y=5_parametric, col sep=comma] {MAE_test_m20_MMNL_n10_L5.csv};
\addplot table [x=MMNL_segments_true, y=10_parametric, col sep=comma] {MAE_test_m20_MMNL_n10_L5.csv};
\end{axis}

\begin{axis}[
name=m50,
at=(m10.below south), anchor=above north west,
ybar,
bar width=\barwid,
width=\wid\textwidth,
height=\hgt\textwidth,
symbolic x coords={1,5,10},
xtick=data,
yticklabel style={/pgf/number format/fixed},
ymin=0,ymax=0.05,
xlabel={$K$},
ylabel={Test set MAE},
legend style={legend pos=outer north east,
			  cells={anchor=west}},
title={$m=50$},
enlarge x limits=0.25,
]
\addplot table [x=MMNL_segments_true, y=1_nonparametric, col sep=comma] {MAE_test_m50_MMNL_n10_L5.csv};
\addplot table [x=MMNL_segments_true, y=2_nonparametric, col sep=comma] {MAE_test_m50_MMNL_n10_L5.csv};
\addplot table [x=MMNL_segments_true, y=inf_nonparametric, col sep=comma] {MAE_test_m50_MMNL_n10_L5.csv};
\addplot table [x=MMNL_segments_true, y=1_parametric, col sep=comma] {MAE_test_m50_MMNL_n10_L5.csv};
\addplot table [x=MMNL_segments_true, y=3_parametric, col sep=comma] {MAE_test_m50_MMNL_n10_L5.csv};
\addplot table [x=MMNL_segments_true, y=5_parametric, col sep=comma] {MAE_test_m50_MMNL_n10_L5.csv};
\addplot table [x=MMNL_segments_true, y=10_parametric, col sep=comma] {MAE_test_m50_MMNL_n10_L5.csv};
\legend{$\ell_1$, $\ell_2$, $\ell_\infty$, $K_{\fit}=1$, $K_{\fit}=3$, $K_{\fit}=5$, $K_{\fit}=10$}
\end{axis}
\end{tikzpicture}
\caption{Test set MAE for $L=5$.}
\label{fig:MAE-par_vs_nonpar}
\end{figure}

comparefig4.tex
\renewcommand{\wid}{0.45}
\renewcommand{\hgt}{0.35}
\renewcommand{\scalefactor}{1}

\begin{figure}[htb]
\centering
\subfigure[Test set MAE]{
\begin{tikzpicture}[scale=\scalefactor]
\begin{axis}[
name=dynamicMAE,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
xticklabels={50,100,200,300,400,500,600,700,800,900,1000,$\infty$},
xticklabel style={font=\footnotesize, rotate=55},
yticklabel style={/pgf/number format/fixed},
ymin=0,ymax=0.029,
xlabel={$\kappa$},
ylabel={Test set MAE},
title={},
enlarge x limits=0.1,
]
\addplot table [x=labels, y=MAE_test, col sep=comma] {dynamic_init2000_MAE_test_MMNL_n10_L5_K5.csv};
\end{axis}
\end{tikzpicture}
\label{fig:dynamic-run}
}%
\hfill
\subfigure[Iterations]{
\begin{tikzpicture}[scale=\scalefactor]
\begin{axis}[
name=dynamicMAE,
width=\wid\textwidth,
height=\hgt\textwidth,
xtick=data,
xticklabels={50,100,200,300,400,500,600,700,800,900,1000,$\infty$},
xticklabel style={font=\footnotesize, rotate=55},
yticklabel style={/pgf/number format/fixed},
ymin=0,ymax=1200,
xlabel={$\kappa$},
ylabel={Iterations},
title={},
enlarge x limits=0.1,
]
\addplot table [x=labels, y=iterations, col sep=comma] {dynamic_init2000_MAE_test_MMNL_n10_L5_K5.csv};
\end{axis}
\end{tikzpicture}
\label{fig:dynamic-increments}
}
\caption{Dynamically updating $p^t$ for $\ell_2$-norm, $L=5$, $K=5$, and $m=20$.}
\label{fig:dynamic}
\end{figure}

Def.tex
\def\crd{\color{red}}
\def\cbl{\color{blue}}
\def\cgr{\color{green}}
\def\cmg{\color{magenta}}
\def\ccy{\color{cyan}}
\def\cye{\color{yellow}}
\definecolor{gold}{rgb}{0.85,0.65,0}

\newcommand{\wrt}{w.r.t.}

\newcommand{\nin}{\noindent}
\newcommand{\nn}{\nonumber \\}
\renewcommand{\theequation}{\mbox{\arabic{equation}}}
\newcommand{\ese}{\end{eqnarray*}}
\newcommand{\bse}{\begin{eqnarray*}}
\def\beq{\begin{equation}}
\def\eeq{\end{equation}}
\def\Bbb#1{{\bf #1}}
\def\eqnum#1{\eqno (#1)}
\def\fnote#1{\footnote}
\newcommand{\epr}{\hfill\hbox{\hskip 4pt \vrule width 5pt height 6pt depth 1.5pt}\vspace{0.0cm}\par}
\def\qed{\ \hfill$\square$\par\smallskip}


\def\ra{\rangle}
\def\la{\langle}
\def\ul{\underline}

\newcommand{\norm}[1]{\left\Vert #1 \right\Vert}

\newcommand{\grad}{\ensuremath{\nabla}}



\def\A{{\mathbb{A}}}
\def\B{{\mathbb{B}}}
\def\C{{\mathbb{C}}}
\def\D{{\mathbb{D}}}
\def\E{{\mathbb{E}}}
\def\F{{\mathbb{F}}}
\def\G{{\mathbb{G}}}
\def\H{{\mathbb{H}}}
\def\I{{\mathbb{I}}}
\def\J{{\mathbb{J}}}
\def\K{{\mathbb{K}}}
\def\L{{\mathbb{L}}}
\def\M{{\mathbb{M}}}
\def\N{{\mathbb{N}}}
\def\O{{\mathbb{O}}}
\def\P{{\mathbb{P}}}
\def\Q{{\mathbb{Q}}}
\def\R{{\mathbb{R}}}
\def\S{{\mathbb{S}}}
\def\T{{\mathbb{T}}}
\def\U{{\mathbb{U}}}
\def\V{{\mathbb{V}}}
\def\W{{\mathbb{W}}}
\def\X{{\mathbb{X}}}
\def\Y{{\mathbb{Y}}}
\def\Z{{\mathbb{Z}}}
\def\Le{{\mathbb{L}}}
\def\Se{{\mathbb{S}}}


\def\bA{{\mathbf{A}}}
\def\bB{{\mathbf{B}}}
\def\bC{{\mathbf{C}}}
\def\bD{{\mathbf{D}}}
\def\bE{{\mathbf{E}}}
\def\bF{{\mathbf{F}}}
\def\bG{{\mathbf{G}}}
\def\bH{{\mathbf{H}}}
\def\bI{{\mathbf{I}}}
\def\bJ{{\mathbf{J}}}
\def\bK{{\mathbf{K}}}
\def\bL{{\mathbf{L}}}
\def\bM{{\mathbf{M}}}
\def\bN{{\mathbf{N}}}
\def\bO{{\mathbf{O}}}
\def\bP{{\mathbf{P}}}
\def\bQ{{\mathbf{Q}}}
\def\bR{{\mathbf{R}}}
\def\bS{{\mathbf{S}}}
\def\bT{{\mathbf{T}}}
\def\bU{{\mathbf{U}}}
\def\bV{{\mathbf{V}}}
\def\bW{{\mathbf{W}}}
\def\bX{{\mathbf{X}}}
\def\bY{{\mathbf{Y}}}
\def\bZ{{\mathbf{Z}}}

\def\cA{\mathcal{A}}
\def\cB{{\cal B}}
\def\cC{\mathcal{C}}
\def\cD{\mathcal{D}}
\def\cE{{\cal E}}
\def\cF{{\cal F}}
\def\cG{{\cal G}}
\def\cH{{\cal H}}
\def\cI{{\cal I}}
\def\cJ{{\cal J}}
\def\cK{{\cal K}}
\def\cL{{\cal L}}
\def\cM{{\cal M}}
\def\cN{{\cal N}}
\def\cO{{\cal O}}
\def\cP{{\cal P}}
\def\cQ{{\cal Q}}
\def\cR{{\cal R}}
\def\cS{\mathcal{S}}
\def\cT{{\cal T}}
\def\cU{{\cal U}}
\def\cV{\mathcal{V}}
\def\cW{{\cal W}}
\def\cX{{\cal X}}
\def\cY{{\cal Y}}
\def\cZ{\mathcal{Z}}

\newcommand{\bbA}{\mathbb{A}}
\newcommand{\bbB}{\mathbb{B}}
\newcommand{\bbC}{\mathbb{C}}
\newcommand{\bbD}{\mathbb{D}}
\newcommand{\bbE}{\mathbb{E}}
\newcommand{\bbF}{\mathbb{F}}
\newcommand{\bbG}{\mathbb{G}}
\newcommand{\bbH}{\mathbb{H}}
\newcommand{\bbI}{\mathbb{I}}
\newcommand{\bbJ}{\mathbb{J}}
\newcommand{\bbK}{\mathbb{K}}
\newcommand{\bbL}{\mathbb{L}}
\newcommand{\bbM}{\mathbb{M}}
\newcommand{\bbN}{\mathbb{N}}
\newcommand{\bbO}{\mathbb{O}}
\newcommand{\bbP}{\mathbb{P}}
\newcommand{\bbQ}{\mathbb{Q}}
\newcommand{\bbR}{\mathbb{R}}
\newcommand{\bbS}{\mathbb{S}}
\newcommand{\bbT}{\mathbb{T}}
\newcommand{\bbU}{\mathbb{U}}
\newcommand{\bbV}{\mathbb{V}}
\newcommand{\bbW}{\mathbb{W}}
\newcommand{\bbX}{\mathbb{X}}
\newcommand{\bbY}{\mathbb{Y}}
\newcommand{\bbZ}{\mathbb{Z}}

\newcommand{\CA}{\mathcal{A}}
\newcommand{\CB}{\mathcal{B}}
\newcommand{\CC}{\mathcal{C}}
\newcommand{\CD}{\mathcal{D}}
\newcommand{\CE}{\mathcal{E}}
\newcommand{\CF}{\mathcal{F}}
\newcommand{\CG}{\mathcal{G}}
\newcommand{\CH}{\mathcal{H}}
\newcommand{\CI}{\mathcal{I}}
\newcommand{\CJ}{\mathcal{J}}
\newcommand{\CK}{\mathcal{K}}
\newcommand{\CL}{\mathcal{L}}
\newcommand{\CM}{\mathcal{M}}
\newcommand{\CN}{\mathcal{N}}
\newcommand{\CO}{\mathcal{O}}
\newcommand{\CP}{\mathcal{P}}
\newcommand{\CQ}{\mathcal{Q}}
\newcommand{\CR}{\mathcal{R}}
\newcommand{\CS}{\mathcal{S}}
\newcommand{\CT}{\mathcal{T}}
\newcommand{\CU}{\mathcal{U}}
\newcommand{\CV}{\mathcal{V}}
\newcommand{\CW}{\mathcal{W}}
\newcommand{\CX}{\mathcal{X}}
\newcommand{\CY}{\mathcal{Y}}
\newcommand{\CZ}{\mathcal{Z}}



\def\mones{{\mathbf{1}}}

\newcommand{\ihat}{\hat{\imath}}
\newcommand{\jhat}{\hat{\jmath}}

\def\Ic{{\overline{I}}}

\def\e{\mathop{{\rm e}}}
\def\nrm{\mathop{{\parallel}}}
\def\vrt{\mathop{\,\mid}\,}
\def\eqbd{\colon =}

\DeclareMathOperator{\Opt}{Opt}
\def\Argmin{\mathop{\rm Argmin}}
\def\Argmax{\mathop{\rm Argmax}}

\DeclareMathOperator*{\argmin}{arg\,min}
\DeclareMathOperator*{\argmax}{arg\,max}

\DeclareMathOperator{\SadVal}{SV}
\DeclareMathOperator{\SV}{SV}
\DeclareMathOperator{\sad}{sad}
\newcommand{\VI}{VI}
\DeclareMathOperator{\Proj}{Proj}
\DeclareMathOperator{\Prox}{Prox}
\DeclareMathOperator{\Entr}{Entr}
\def\epsilonsad{\epsilon_{\sad}}
\def\epsilonvi{\epsilon_{\hbox{\scriptsize VI}}}
\def\epsilonN{\epsilon_{\hbox{\scriptsize Nash}}}
\def\Sad{{\mathop{\rm Sad}}}
\def\Residual{\mathop{{\rm Residual}}}
\def\DualityGap{{\hbox{\rm DualityGap}}}
\newcommand{\dgf}{d.g.f.}



\DeclareMathOperator{\sign}{sign}
\DeclareMathOperator{\sgn}{sgn}
\DeclareMathOperator{\Supp}{supp}
\DeclareMathOperator{\supp}{supp}
\DeclareMathOperator{\Card}{Card}
\DeclareMathOperator{\card}{card}
\DeclareMathOperator{\Size}{Size}
\DeclareMathOperator{\size}{Size}
\DeclareMathOperator{\Sz}{Size}
\DeclareMathOperator{\Vol}{Vol}
\DeclareMathOperator{\vol}{Vol}
\DeclareMathOperator{\Diam}{Diam}
\DeclareMathOperator{\AvDiam}{AvDiam}
\def\const{{\hbox{\rm const}}}
\def\data{\hbox{\rm data}}
\def\nnz{{\hbox{\rm nnz}}}

\DeclareMathOperator{\rank}{rank}
\DeclareMathOperator{\Diag}{Diag}
\DeclareMathOperator{\diag}{diag}
\DeclareMathOperator{\Trace}{Trace}
\DeclareMathOperator{\Tr}{Tr}
\DeclareMathOperator{\Det}{Det}
\DeclareMathOperator{\Col}{Col}
\DeclareMathOperator{\Row}{Row}
\DeclareMathOperator{\Null}{Null}
\DeclareMathOperator{\Image}{Im}
\DeclareMathOperator{\Imm}{Im}
\DeclareMathOperator{\Range}{Range}
\def\Arrow{{\hbox{\rm Arrow}}}
\def\Fr{{\hbox{\scriptsize Fr}}}
\def\Fro{\mathop{\hbox{\rm\scriptsize Fro}}}
\def\Spec{\mathop{\hbox{\rm\scriptsize Spec}}}


\def\Prob{{\hbox{\rm Prob}}}
\def\Cov{{\hbox{\rm Cov}}}
\def\uniform{\mathop{\rm uniform}}
\DeclareMathOperator{\Uniform}{Uniform}
\DeclareMathOperator{\Geometric}{Geometric}
\DeclareMathOperator{\Poisson}{Poisson}
\DeclareMathOperator{\Var}{Var}
\DeclareMathOperator{\CVaR}{CVaR}
\DeclareMathOperator{\KL}{KL}
\DeclareMathOperator{\fit}{fit}


\def\dim{\mathop{{\rm dim}\,}}
\DeclareMathOperator{\codim}{codim}
\DeclareMathOperator{\Lin}{Lin}
\DeclareMathOperator{\lin}{lin}
\DeclareMathOperator{\Aff}{Aff}
\DeclareMathOperator{\aff}{aff}
\DeclareMathOperator{\Span}{Span}
\DeclareMathOperator{\intt}{int}
\DeclareMathOperator{\myint}{int}
\DeclareMathOperator{\inter}{int}
\DeclareMathOperator{\relint}{rel int}
\DeclareMathOperator{\rint}{rel int}
\DeclareMathOperator{\ri}{rel int}
\DeclareMathOperator{\partialri}{rel int}
\DeclareMathOperator{\cl}{cl}
\DeclareMathOperator{\bd}{bd}
\DeclareMathOperator{\rbd}{rbd}

\DeclareMathOperator{\Dom}{Dom}
\DeclareMathOperator{\dom}{dom}
\DeclareMathOperator{\Epi}{Epi}
\DeclareMathOperator{\epi}{epi}
\DeclareMathOperator{\Ker}{Ker}

\DeclareMathOperator{\Conv}{Conv}
\DeclareMathOperator{\conv}{conv}
\DeclareMathOperator{\co}{conv}
\DeclareMathOperator{\Cone}{Cone}
\DeclareMathOperator{\cone}{cone}
\DeclareMathOperator{\clconv}{\overline{conv}}
\DeclareMathOperator{\clcone}{\overline{cone}}
\DeclareMathOperator{\ccvh}{\overline{conv}}
\DeclareMathOperator{\ccnh}{\overline{cone}}
\def\Ext{{\mathop{\rm Ext}}}
\def\Extr{{\mathop{\rm Extr}}}
\DeclareMathOperator{\Rec}{Rec}
\DeclareMathOperator{\rec}{Rec}
\DeclareMathOperator{\apex}{apex}


\def\dist{\mathop{{\rm dist}}}
\def\log{\mathop{{\rm log}}}
\def\asin{{\hbox{\rm asin}}}
\def\Polar{{\mathop{{\rm Polar\,}}}}
\def\ch{{\mathop{\rm ch\,}}}
\def\arc{{\mathop{\rm arc\,}}}
\def\EllIn{{\mathop{\rm EllIn}}}
\def\GDA{{\rm \ ArD\ }}
\def\GDS{{\rm \ StD\ }}
\def\lev{{\mathop{\rm lev}}}


\def\rdiamondsuit{{\color{red}$\diamondsuit$}}
\def\rclubsuit{{\color{red}$\clubsuit$}\ }
\def\rspadesuit{{\color{red}$\spadesuit$}\ }
\def\rbullet{{\color{red}$\bullet$}\ }
\def\rheartsuit{{\color{red}$\heartsuit$}}
\def\rArrow{{\color{red}$\Rightarrow$}\ }
\def\blacksquare{\hbox{\vrule width 4pt height 4pt depth 0pt}}
\def\square{\hbox{\vrule\vbox{\hrule\phantom{o}\hrule}\vrule}}


\def\nr{\par \noindent}
\def\pn{\par\noindent}
\def\qr{\qquad\cr}
\def\fig{\epsfxsize=220pt\epsfysize=140pt\epsffile}
\def\mypictsmall{\epsfxsize=120pt\epsfysize=120pt\epsffile}
\def\mypictvsmall{\epsfxsize=90pt\epsfysize=90pt\epsffile}

\def\huh{\hbox{\vrule width 2pt height 8pt depth 2pt}}



\def\CE{\mathop{{\hbox{\sl CircEll}}}}
\def\CS{\mathop{{\hbox{\sl CircSimp}}}}
\def\IE{\mathop{{\hbox{\sl InscrEll}}}}
\def\Compl{\mathop{{\rm Compl}}}
\def\Accur{\mathop{{\rm Accur}}}
\def\mes{\mathop{{\rm mes}}}
\def\EllOut{\mathop{{\rm EllOut}}}
\def\Ell{\mathop{{\sl Ell}}}
\def\InnerRad{\mathop{{\rm InnerRad}}}
\def\Mom{\mathop{\rm Mom}}
\def\Loss{{\hbox{\rm Loss}}}


\def\bStnu{\bS^{\widetilde{\nu}}}
\def\bStnup{\bS^{\widetilde{\nu}}_+}




DynamicLearning.tex
\section{Dynamic Estimation of a Non-Parametric Choice Model}\label{sec:learning-choice-model-MD}


Our general approach to learning a non-parametric choice model is based on inferring $\lambda$ by minimizing a distance measure $D(\cdot,\cdot)$ from the theoretical probabilities $A\lambda = x \in X$ to the empirical observations $p$:
\begin{equation}\label{eqn:general-distance-learn}
\min_{\lambda} \left\{ D(A\lambda,p):~ \lambda \in \Delta_{n!} \right\} = \min_x \left\{ D(x,p):~ x \in X \right\}.
\end{equation}
We give a brief overview of three closely related approaches to learn the non-parametric choice model from data \citep{FariasJS2013,vanRyzinVulcano2015,BertsimasMisic2015} in \ref{sec:learning-choice-model}. There, we also demonstrate that these seemingly disparate models used in the literature
are special instantiations of this more general distance minimization framework.

In this paper, we assume that $D(\cdot,p)$ is convex and continuous on its domain, and that $D(x,x)=0$; this is the case for the specific instantiations of $D$ in \citet{FariasJS2013,vanRyzinVulcano2015,BertsimasMisic2015} as well.

In the static setup, for a specified accuracy level $\epsilon>0$, our goal is to obtain an (additive error) $\epsilon$-approximate solution to \eqref{eqn:general-distance-learn} within a reasonable number of iterations. The main challenge in solving \eqref{eqn:general-distance-learn} is that in general the set $X$ only admits a high-dimensional representation as $X = \{x = A\lambda :~ \lambda \in \Delta_{n!}\}$. Thus, projection-based or interior point methods will not work for $X$. On the other hand, \emph{linear} optimization over $X$, while non-trivial, is a manageable problem. We discuss this in Section \ref{sec:comb-subproblem}. Furthermore, in order to compute choice probabilities for unseen subsets, when given a candidate solution $x$, we also need to identify the rankings $\sigma$ that generated $x$, i.e., we want to know $\lambda$ such that $x=A\lambda$. The solution frameworks that we present aim to circumvent the challenging representation of $X$, while providing a candidate solution in terms of the vertices $a(\sigma)$.

In the dynamic variant of the problem, instead of fixed data $p$, we now have changing data $p^t$ that converges to a limit $p$ as we have more observations, i.e., $t\to\infty$.
Under standard statistical assumptions, the estimates $p^t$ obtained through \eqref{eqn:pq_{ij}} from a growing set of observations converge to the true distribution vector $p = \{\bbP[i, \mid \CA_j]\}_{i \in \CA_j, j \in [m]}$ (almost surely).
In this setup we would still like to solve \eqref{eqn:general-distance-learn}, but we are only given access to the sequence $\{p^t\}_{t \geq 1}$. In particular, we are interested in finding a solution $x^T$ after $T$ iterations (or time steps) along with a bound on its error, i.e.,  $D\left( x^T, p \right) - \min_{x \in X} D(x,p)$, that converges to zero as the number of iterations $T$ goes to infinity. 
This is now exactly a joint estimation and optimization (JEO) problem studied in \citet{AhmadiShanbhag2014,HoNguyenKK2016}, for which methods exist to optimize $f(y,p)$ using a sequence $p^t \to p$ and employing efficient updates at each step $t$. Our two solution frameworks for \eqref{eqn:general-distance-learn} can be viewed as two new methods to solve the JEO problem designed specifically to handle the high-dimensionality issues in  problems such as non-parametric choice model estimation.



\subsection{Primal Approach via the Frank-Wolfe Algorithm}\label{sec:FW}

We first discuss solving \eqref{eqn:general-distance-learn} via the Frank-Wolfe (F-W) algorithm under a finite curvature constant assumption. In the context of the static setup, \citet{Jaggi2013} derives the standard F-W error bound of $O(1/T)$ after $T$ iterations. However, this is explicitly based on the finiteness of an important quantity referred to as the \emph{curvature constant} (see \eqref{eqn:curvature-constant}), and we show that when $D$ is a $\ell_q$-norm or weighted KL-divergence, the curvature constant is infinite; so the usual F-W bounds cannot apply. Nevertheless, when the distance measure $D$ is set to be the  square of an $\ell_q$ norm, $D$ satisfies $1$-smoothness, and in such a case we can  guarantee a rate of $O(1/\sqrt{T})$ convergence for $\ell_q$-norm minimization. Moreover, in the dynamic data setting, we establish that a natural variant of the F-W algorithm also enjoys convergence guarantees under an additional minimum data convergence rate assumption. 

In addition to our general convexity assumption on $D(\cdot,p)$, our primal approach is based on the following regularity condition, which can be viewed as a generalized triangle inequality.
\begin{assumption}\label{ass:primal-distance-regular}
For some increasing continuous function $g_D:[0,\infty) \to [0,\infty)$ such that $g(0)=0$, $D(x,p) - D(x,p') \leq g_D(D(p',p))$ for all $x,p,p'$.
\end{assumption}
Following \citet{Jaggi2013}, we define the curvature constant critical in analysing F-W methods and assume it is finite.
\begin{assumption}\label{ass:distance-smoothness}
The \emph{curvature constant} of $D$ defined below is finite and uniformly bounded in $p^t$: 
\begin{equation}\label{eqn:curvature-constant}
C_{D,t} := \sup_{\substack{x,s \in X\\ \alpha \in [0,1]}} \frac{1}{\alpha^2} \left( D((1-\alpha)x+\alpha s,p^t) - D(x,p^t) - \alpha \la s-x, \grad_x D(x,p^t) \ra \right)
\end{equation}
and $C_{D,t} \leq C_D < \infty$ for all $t \geq 1$. 
\end{assumption}

The variant of F-W algorithm for the dynamic setup is stated in Algorithm \ref{alg:FW}. The key difference of Algorithm \ref{alg:FW} as opposed to the standard F-W algorithm is that each step $t$, Algorithm \ref{alg:FW} works with a gradient of a dynamically changing $D(\cdot,p^t)$ instead of a fixed $D(\cdot,p)$. This setup is similar to the online F-W algorithms; see e.g.,  \citet{HazanKale2012}. Note that \citet{HazanKale2012} and other online algorithms are usually concerned with obtaining regret bounds. On the other hand, our analysis of Algorithm \ref{alg:FW} in the dynamic setup is not based on such regret bounds.
\begin{algorithm}[htb]
\caption{Frank-Wolfe (F-W) algorithm for solving \eqref{eqn:general-distance-learn}.}\label{alg:FW}
\begin{algorithmic}
\STATE {\bf Input:~} time horizon $T$, initial point $x^1 \in X$, step sizes $\{\gamma^t\}_{t=1}^T$, $\gamma_t \in [0,1]$.
\FOR{$t=1,\ldots,T-1$}
	\STATE $z^t = \argmin_{z \in X} \la \grad_x D(x^t,p^t), z \ra$.
	\STATE $x^{t+1} = (1-\gamma^t) x^t + \gamma^t z^t$.
\ENDFOR
\STATE {\bf Output:~} solution $x^T \in X$.
\end{algorithmic}
\end{algorithm}

Since $X$ is a polytope, and $z^t$ is the result of a linear optimization over $X$, $z^t$ will be a vertex, i.e., a column $a(\sigma^t)$ for some ranking $\sigma^t \in S_n$. Thus, the output $x^T$ corresponds to a distribution $\lambda^T \in \Delta_{n!}$, with sparsity at most $T$. We discuss the linear optimization subproblem arising in Algorithm \ref{alg:FW} for the computation of $z^t$ further in Section \ref{sec:comb-subproblem}.

We next derive the convergence guarantee of Algorithm \ref{alg:FW} in the dynamic setup under our assumptions; see \ref{sec:proofs} for the proofs in this section.
\begin{theorem}\label{thm:FW-dynamic-convergence}
Suppose that Assumptions \ref{ass:primal-distance-regular} and \ref{ass:distance-smoothness} hold, and that Algorithm \ref{alg:FW} is ran for $T\geq 4$ iterations to find a point $x^T \in X$ with step sizes $\gamma_t = 2/(t+1)$. Then for any $x \in X$,
\begin{align*}
&D(x^T,p) - D(x,p)\\
&\leq \frac{4C_D}{T} + \frac{1}{(T-1)T} \sum_{t=1}^{T-3} t(t+1) g_D(D(p^t,p)) + \frac{1}{(T-1)T} \sum_{t=1}^{T-3} \left( (t-1)t + \frac{4t}{t+1}\right) g_D(D(p,p^t))\\
&\quad + \frac{T-2}{T} \left( \frac{T-3}{T-1}  + \frac{4}{(T-1)^2} \right) g_D(D(p,p^{T-2})) + \left(\frac{T-2}{T} + \frac{4}{T^2}\right) g_D(D(p,p^{T-1})).
\end{align*}
\end{theorem}
Notice in Theorem \ref{thm:FW-dynamic-convergence} that the error bound has the $4C_D/T$ rate along with four error terms which are  accumulated from using approximate data $p^t \approx p$. Indeed, if $p^t = p$ for all $t$, i.e., the static case, these additional error terms disappear, and we are left with only the $4C_D/T$ term. In order to get convergence in the dynamic case, these error terms must go to $0$ as $T \to \infty$. To ensure this, we must require that $D(p^t,p)$ converges sufficiently fast to guarantee that $\sum_{t=1}^{T-3} t(t+1) g_D(D(p^t,p)) = o(T^2)$. This is a major handicap of the primal F-W based approach. Our dual-based approach outlined in Section \ref{sec:saddle-point} avoids such a requirement, and instead only needs $D(p^t,p), D(p,p^t) \to 0$.

Note that in this approach, the overall number of iterations $T$ upper-bounds the sparsity of the learned choice model $\lambda$ as well. As a result, our analysis in essence exposes an explicit trade-off between the sparsity of the non-parametric choice model $\lambda$ and its estimation accuracy $\epsilon$.

There are a number of different instantiations of the distance measure $D$ suggested in the literature. Two such choices are the weighted KL-divergence (which stems from the maximum likelihood) and $\ell_1$-norm suggested in \citet{vanRyzinVulcano2015} and \citet{BertsimasMisic2015} respectively. See \ref{sec:learning-choice-model} for precise details on how the approaches in \citet{vanRyzinVulcano2015,BertsimasMisic2015} correspond to particular choices of $D$ in our framework. Unfortunately, for both of these choices of $D$, in Proposition~\ref{prop:infinite-curvature constant} we demonstrate that the curvature constant \eqref{eqn:curvature-constant}, which is critical in terms of ensuring the convergence of F-W algorithm, is infinite. Thus, this rules out the use of the classical F-W algorithm in the static case and the application of Theorem \ref{thm:FW-dynamic-convergence} in the dynamic case when $D$ is selected to be either an $\ell_q$ norm or the weighted KL-divergence.
\begin{proposition}\label{prop:infinite-curvature constant}
Suppose $n>2$. 
For any $q \in [1,\infty]$, the function $D(x,p) = \|x-p\|_q$ has infinite curvature constant \eqref{eqn:curvature-constant} for any $p \in X$. Furthermore, when $D(x,p)=\sum_{j \in [m]} w_j \KL(p_j,x_j)$ for any positive weights $w_j$, and $p \in X$ such that $p_{ij} > 0$ for all $i \in \CA_j$, the curvature constant is infinite.
\end{proposition}

As opposed to Proposition~\ref{prop:infinite-curvature constant}, for $D(x,p) = \|x-p\|^2$ for some norm $\|\cdot\|$, it is possible to bound the associated curvature constant and Assumption~\ref{ass:primal-distance-regular} holds; we discuss this below.

\begin{remark}\label{rem:norm-FW}
Suppose we define $D(x,p) = \|x-p\|$ for some norm $\|\cdot\|$. Proposition \ref{prop:infinite-curvature constant} states that the usual F-W bounds cannot apply if we use $D(x,p) = \|x-p\|_q$. However, note that the function $\frac{1}{2} \|\cdot\|_*^2$ is (trivially) 1-strongly convex with respect to the dual norm $\|\cdot\|_*$. Moreover, the convex conjugate of $\frac{1}{2} \|\cdot\|_*^2$ is $\frac{1}{2} \|\cdot\|^2$, and hence the standard convex analysis results imply that it is $1$-smooth with respect to $\|\cdot\|$. Thus, as stated in \citet{Jaggi2013}, for $D(x,p) = \frac{1}{2} \|x-p\|^2$, we can bound the associated curvature constant $C_D$ by the diameter $R = \max_{x,x' \in X} \|x-x'\|$ of $X$, i.e., $C_D \leq R^2$. Since $X$ is a polytope, and hence compact, $R < \infty$. We can then apply Theorem \ref{thm:FW-dynamic-convergence} to ensure that $\|x^T-p\|^2 \leq \min_{x \in X} \|x-p\|^2 + 8R^2/T$ in the static case, and by taking the square root, we have $\|x^T-p\| - \min_{x \in X} \|x-p\| \leq 2\sqrt{2} R/\sqrt{T}$. In the dynamic case with changing data $p^t$ at each time step $t$, the norm function satisfies $\frac{1}{2} \|x-p\|^2 - \frac{1}{2} \|x-p^t\|^2 = \frac{1}{2} (\|x-p\|-\|x-p^t\|)(\|x-p\|+\|x-p^t\|) \leq R \|p-p^t\| = R \sqrt{2 \cdot \frac{1}{2}\|p-p^t\|^2}$ due to the reverse triangle inequality, so Assumption \ref{ass:primal-distance-regular} holds with $g_D(w) = R \sqrt{2w}$ and convergence is guaranteed when $p^t \to p$ sufficiently fast.
\halmos
\end{remark}



\subsection{Dual Approach via Regret Minimization}\label{sec:saddle-point}

A major handicap of the primal approach based on F-W algorithm is that in the dynamic case the overall convergence rate depends heavily on the rate of convergence for $p^t \to p$. In this section, we explore an alternative dual approach  under Assumption \ref{ass:dual-distance-regular} that removes such a dependence completely. For the dual approach, we define the convex conjugate
\[D_*(y,p) := \sup_z \left\{ \la y,z \ra - D(z,p) :~ z \in \bbR^N \right\}\]
of $D(\cdot,p)$, and assume the following regularity conditions on $D$.
\begin{assumption}\label{ass:dual-distance-regular}
All subgradients $\|\grad_x D(x,p)\| \leq G$ of $D(\cdot,p)$ are bounded for any $p$ and some norm $\|\cdot\|$. Furthermore, for some $0 \leq L_D, L_{D,*} < \infty$, $\min_{x \in X} D(x,p) - \min_{x \in X} D(x,p') \leq L_D D(p,p')$ and $D_*(y,p) - D_*(y,p') \leq L_{D,*} D(p',p)$ for all $y,p,p'$.
\end{assumption}
When $D(x,p) = \|x-p\|$ for any norm, we discuss in Remark~\ref{rem:norm-assumptions-example} that Assumption~\ref{ass:dual-distance-regular} is satisfied.

Given that $D(\cdot,p)$ is convex and continuous on its domain, we can write $D(x,p)$ via its convex conjugate $D_*$:
\begin{align*}
D(x,p) &= \sup_y \left\{ \la x,y \ra - D_*(y,p) :~ y \in \bbR^N \right\}.
\end{align*}
Note that $D_*(y,p)$ is convex in $y$. By Assumption \ref{ass:dual-distance-regular}, $\|\grad_x D(x,p)\| \leq G$, so denoting $Y := \{y :~ \|y\| \leq G\}$ we can restrict the domain of $y$ accordingly:
\[ D(x,p) = \sup_{y \in Y} \left\{ \la x,y \ra - D_*(y,p) :~ y \in \bbR^N \right\}. \]
Based on these definitions, the problem \eqref{eqn:general-distance-learn} now admits a natural saddle point representation:
\begin{equation}\label{eqn:norm-min-SP}
\SV(p) := \min_{x} \left\{ D(x,p):~ x \in X \right\} = \min_{x \in X} \sup_{y \in Y} \left\{ \la x,y \ra - D_*(y,p) \right\}.
\end{equation}
Ideally, we would solve \eqref{eqn:norm-min-SP} with a saddle point algorithm such as Mirror Prox \citep{Nemirovski2004}, which leads to a convergence rate of $O\left( 1/T \right)$ after $T$ iterations. However, each iteration of Mirror Prox, or other specialized saddle point algorithms,  involves projection operations, which can be difficult for the polytope $X$. Indeed, in such schemes, one needs to update every component of $\lambda$, the underlying probability distribution over rankings $\sigma \in S_n$, which has $n!$ entries and is intractable.

To remedy this, we directly solve the dual problem
\begin{equation}\label{eqn:norm-min-convex-form}
\SV(p) = \sup_{y \in Y} f(y,p),
~~~\text{where}~~~
f(y,p) := \min_{\lambda \in \Delta_{n!}} \left\{ \la A\lambda,y \ra - D_*(y,p) \right\} 
= \min_{x\in X} \left\{ \la x,y \ra - D_*(y,p) \right\} 
. 
\end{equation} 
Notice that $f(y,p)$ is concave in $y$ and has supergradients $\grad_y f(y,p) = x \!-\! \grad_y D_*(y,p)$, where $x \in \argmin_{x' \in X} \! \la x',y \ra$. Thus, $\SV(p)$ is simply maximizing a concave function over a convex norm ball $Y$, which is much more manageable than $X$.

In the static case, we can recover a primal solution $x \in X$ for the original problem \eqref{eqn:general-distance-learn} via solving the dual problem \eqref{eqn:norm-min-convex-form} by using an approach from \citet{NedicOzdaglar2009}. 
We note that similar ideas were applied to online stochastic programming \citep{AgrawalDevanur2015} and in an algorithmic approach to Approximate Caratheodory Theorem \citep{MirrokniLemeVladuWong2017}. 
As opposed to these approaches in the static setup, in our dynamic setting we only have access to approximate data $p^t$, and hence we only see an approximate sequence of functions $f(y,p^t)$. 

We now discuss how to extend these methods to recover a primal solution $x \in X$ for \eqref{eqn:general-distance-learn} via approximate dual functions $f(y,p^t) \approx f(y,p)$.



The key idea is as follows. For a given sequence of dual points $\{y^t\}_{t=1}^T$, we define concave functions $f^t$ based on the current $y^t$ and $p^t$ by
\begin{equation}\label{eqn:f^t-def}
f^t(y) := \la z^t, y \ra - D_*(y,p^t), \quad \text{where} \quad z^t := \argmin_{z \in X} \la z, y^t\ra.
\end{equation}
Notice in \eqref{eqn:f^t-def} that our functions $f^t(\cdot)$ are defined based on a sequence of primal points $\{z^t\}_{t=1}^T$. From these primal points $\{z^t\}_{t=1}^T$, we build a candidate primal solution $x^T := \frac{1}{T} \sum_{t=1}^T z^t$ and bound its optimality gap. Using these ideas, we derive the following result on the convergence guarantee of this dual approach; see  \ref{sec:proofs} for the corresponding proof.

\begin{theorem}\label{thm:online-choice-model-learning}
Suppose that Assumption \ref{ass:dual-distance-regular} holds. Given a sequence $\{y^t\}_{t=1}^T$, let $\{z^t\}_{t=1}^T$ be generated according to \eqref{eqn:f^t-def}, and $x^T := \frac{1}{T} \sum_{t=1}^T z^t$. Then
\begin{align*}
D\left( x^T, p \right) - \min_{x \in X} D(x,p) &\leq \max_{y \in Y} \frac{1}{T} \sum_{t=1}^{T} f^t(y) - \frac{1}{T} \sum_{t=1}^T f^t(y^t) + \frac{L_D + L_{D,*}}{T} \sum_{t=1}^T D(p^t,p).
\end{align*}
\end{theorem}

Theorem~\ref{thm:online-choice-model-learning} decomposes the optimality gap of our candidate point $x^T$ into two error terms. The latter term $\frac{2}{T} \sum_{t=1}^T D(p^t,p)$ is due to the accumulation of errors from the use of approximate data $p^t \approx p$, and they disappear in the static case when $p^t = p$. In contrast to Theorem \ref{thm:FW-dynamic-convergence}, this term converges to $0$ whenever $D(p^t,p) \to 0$, \emph{regardless} of the rate of convergence of $p^t \to p$. The first term is a regret term for the sequence $\{y^t\}_{t=1}^T$ on the functions $\{f^t\}_{t=1}^T$.
Since $f^t(y) \geq f(y,p^t)$ and $f^t(y^t) = f(y^t,p^t)$, we can bound the regret on the (approximate) dual objectives $\{f(y,p^t)\}_{t=1}^T$ as follows:
\[ \max_{y \in Y} \frac{1}{T} \sum_{t=1}^{T} f(y,p^t) - \frac{1}{T} \sum_{t=1}^T f(y^t,p^t) \leq \max_{y \in Y} \frac{1}{T} \sum_{t=1}^{T} f^t(y) - \frac{1}{T} \sum_{t=1}^T f^t(y^t). \]
Thus, by choosing $\{y^t\}_{t=1}^T$ to have small regret on $\{f^t\}_{t=1}^T$, we are also near-optimal on the dual objectives. Online convex optimization (OCO) gives us several methods for minimizing regret. We state a typical bound below and defer a more complete overview of this and the online Mirror Descent algorithm which achieves this bound to \ref{sec:prox-setup}.
\begin{theorem}
\label{thm:OCO-regret}
Suppose that $\grad_y D_*(y,p^t)$ is bounded for all $y \in Y$ and $p^t$, $t \geq 1$. Then there exists a method to choose $y^t$ using only $y^{t-1},f^{t-1}$ so as to guarantee
\[ \max_{y \in Y} \frac{1}{T} \sum_{t=1}^{T} f^t(y) - \frac{1}{T} \sum_{t=1}^{T} f^t(y^t) \leq O\left( \frac{1}{\sqrt{T}} \right). \]
\end{theorem}
Similar to the primal approach, the overall number of iterations $T$ upper bounds the sparsity of our learned model in the dual approach as well. We now show that the distance measures $D$ based on general norms satisfy Assumption~\ref{ass:dual-distance-regular} required in the dual approach.

\begin{remark}\label{rem:norm-assumptions-example}
Suppose $D(x,p) = \|x-p\|$ is a norm. Standard convex analysis results imply that $\|\grad_y D(y,p)\|_* \leq 1$, hence we define $Y := \{y : \|y\|_* \leq 1\}$ as the unit dual norm ball. Note that $D_*(y,p) = \la y,p \ra$ when $\|y\|_* \leq 1$ and $\infty$ otherwise. For any $y\in Y$, we have $\grad_y D_*(y,p^t) = p^t$ is bounded because the data vectors $p^t$ are bounded. Thus, we can apply Theorem \ref{thm:OCO-regret}. Furthermore, $D_*(y,p^t) - D_*(y,p) = \la y, p^t - p \ra \leq \|p^t - p\|$ for any $y \in Y$. Then from the triangle inequality $\|x-p^t\|\leq \|x-p\| + \|p-p^t\|$, we have $\min_{x \in X} D(x,p^t) - \min_{x \in X} D(x,p) \leq \|p^t-p\|$. Thus Assumption \ref{ass:dual-distance-regular} is satisfied with $L_D = L_{D,*} = 1$. This shows that when $D$ is defined as any norm (in particular the $\ell_1$-norm from \citet{BertsimasMisic2015}), we can estimate a non-parametric choice model from a continuously updated sequence of data $p^t \to p$ by solving a regret minimization problem.
\epr
\end{remark}

Note that, similar to the F-W algorithm, in our dual approach governed by the relations \eqref{eqn:f^t-def}, each $z^t$ is obtained by solving a \emph{linear} optimization problem over $X$. Hence, $z^t$ will be a vertex, i.e., a column $a(\sigma^t)$ for some ranking $\sigma^t \in S_n$, so $x^T$ corresponds to a distribution $\lambda^T$ with sparsity at most $T$. We discuss the optimization subproblem in \eqref{eqn:f^t-def} to compute $z^t$ in Section \ref{sec:comb-subproblem}. 

\subsection{Combinatorial Subproblem}\label{sec:comb-subproblem}

In both approaches of Sections~\ref{sec:FW}~and~\ref{sec:saddle-point}, we must solve a linear optimization problem over $X$. Since $X$ is a polytope with vertices $a(\sigma)$, we have, for a cost vector $c$,
\[ z^* = a(\sigma^*), \quad \text{where} \quad \sigma^* = \argmin_{\sigma \in S_n} \la a(\sigma), c \ra, \]
Thus, in each iteration of the primal and the dual approaches, we must solve the following combinatorial optimization problem over rankings:
\begin{equation}\label{eqn:subgradient-subproblem}
\min_{\sigma} \left\{ \sum_{j \in [m]} \sum_{i \in \cA_j} y_{ij}^t a_{ij}(\sigma) :~ \sigma \in S_n \right\}.
\end{equation}
On the one hand, this problem is NP-hard, since it is a generalization of the linear ordering problem and the maximum weighted independent set problem, see e.g., \cite[Proposition 3]{vanRyzinVulcano2015}. This is the main drawback of our approaches. However, we note that the exact same combinatorial problem must be solved in all other approaches of learning a non-parametric choice model (see \ref{sec:learning-choice-model} and in particular the equations \eqref{eqn:revenue-prediction-robust-constraint}, \eqref{eqn:MLE-column-generating}, \eqref{eqn:norm-l1-min-column-generating}). On the other hand, while we cannot avoid the NP-hardness in learning a non-parametric choice model from data, we note that \eqref{eqn:subgradient-subproblem} can be formulated as a (relatively) compact integer program with $O(n^2)$ variables and $O(n^3)$ constraints, and also it can be handled efficiently by off-the-shelf integer programming solvers. Furthermore, \citet{JagabathulaRusmevichientong2016} prove polynomial-time solvability of \eqref{eqn:subgradient-subproblem} under a number of assumptions on the subsets $\CA_1,\ldots,\CA_m$.
Thus, by employing our suggested first order approaches in Sections~\ref{sec:FW}~and~\ref{sec:saddle-point}, we avoid the problem of having to deal with the high-dimensional representation of $X$ (which can involve $n!$ variables), and instead must solve a relatively compact integer program at each iteration.

Finally, we remark that our dynamic primal and dual approaches can apply to arbitrary domains $X$ besides the ones specified for learning non-parametric choice models in \citet{FariasJS2013,vanRyzinVulcano2015,BertsimasMisic2015}. This flexibility can be attractive in utilizing additional a priori structural information on the choice model $\lambda$.

ExistingApproaches.tex
\section{Existing Approaches to Non-Parametric Choice Estimation}\label{sec:learning-choice-model}

In this section, we examine the existing approaches to learn the non-parametric choice model, i.e., infer an appropriate probability vector $\lambda$ using the data collected via the process outlined in Section~\ref{sec:model-data}, and demonstrate how they are particular instantiations of our general model. 

\subsection{Revenue Prediction Approach}\label{sec:rev-predict}

Let $r_i$ be the revenue of item $i\in[n]$. Then the expected revenue of an assortment $\cA \subset [n]$ under distribution $\lambda$ is $\sum_{i \in \cA} r_i \bbP_{\lambda}[i \mid \cA]$.
\citet{FariasJS2013} seek to find the worst-case expected revenue from a distribution $\lambda$ consistent with the given data in the sense that the theoretical probabilities $\bbP_\lambda[i \mid \cA_j] = \la a_{ij}, \lambda \ra$ are precisely consistent with their empirical estimates $p_{ij}$. Since the probabilities $\bbP_{\lambda}[i \mid \cA]$ are linear in $\lambda$, this can be formulated as a linear program (LP)
\begin{align*}
\min_{\lambda} \left\{ \sum_{i \in \cA} r_i \bbP_{\lambda}[i \mid \cA] :~ A\lambda = p,~ \lambda \in \Delta_{n!}\right\}.
\end{align*}
We first make a few observations related to this model of \citet{FariasJS2013}. 
In fact, when $\cA = \cA_j$ for some $j \in [m]$, we have $\bbP_{\lambda}[i \mid \cA] = \la a_{ij}, \lambda \ra = p_{ij}$ due to the constraints $A\lambda = p$, hence the objective is constant. Thus the LP becomes a feasibility problem
\begin{align}\label{eqn:revenue-prediction-model}
\text{find} \quad \lambda \in \Delta_{n!} \quad \text{s.t.} \quad  A\lambda = p.
\end{align}
That said, \eqref{eqn:revenue-prediction-model} is still computationally intractable even for moderate values of $n$ because it involves $n!$ variables. Nonetheless, the dual of \eqref{eqn:revenue-prediction-model} admits the following robust LP interpretation:
\begin{align}\label{eqn:revenue-prediction-dual}
\max_{\beta,\nu} \left\{ \la \beta, p\ra - \nu :~ \max_{\sigma \in S_n} \la \beta, a(\sigma)\ra \leq \nu \right\}.
\end{align}
Note that verifying the feasibility of a solution with respect to the robust constraint in \eqref{eqn:revenue-prediction-dual}, i.e., 
\begin{equation}\label{eqn:revenue-prediction-robust-constraint}
\max_{\sigma \in S_n} \la \beta, a(\sigma) \ra= \max_{\sigma} \left\{ \sum_{j \in [m]} \sum_{i \in \cA_j} \beta_{ij} a_{ij}(\sigma) :~ \sigma \in S_n \right\} \leq \nu
\end{equation}
is a combinatorial problem of the exact same form as \eqref{eqn:subgradient-subproblem}.
\citet{FariasJS2013} suggests solving \eqref{eqn:revenue-prediction-dual} either using the constraint sampling technique \citep{CalafioreCampi2005} or by building an approximation to its robust counterpart obtained from approximating the uncertainty sets with an efficiently representable polyhedron.

In fact, \eqref{eqn:revenue-prediction-model} can be seen as choosing $\lambda \in \Delta_{n!}$ to minimize a (very harsh) distance measure:
\begin{equation}\label{eqn:revenue-prediction-distance-min}
\min_{\lambda \in \Delta_{n!}} D(A\lambda,p), \quad D(A\lambda,p) = \begin{cases} 0, & A\lambda = p\\ \infty, & \text{otherwise}. \end{cases}
\end{equation}
In general, and specifically when the observations are noisy, there is no guarantee that there exists $\lambda \in \Delta_{n!}$ to fit the data $p$ exactly, i.e., $A\lambda = p$. To remedy this, \citet{vanRyzinVulcano2015} and \citet{BertsimasMisic2015} examine approaches that use less harsh distance measures $D(\cdot,\cdot)$.

\subsection{Maximum Likelihood Estimation Approach}\label{sec:MLE-approach}

\citet{vanRyzinVulcano2015} propose the following method to learn $\lambda$ via maximum likelihood estimation (MLE). We next describe their method and provide an alternative interpretation of their approach as the minimization of a particular distance measure, namely  Kullback-Leibler (KL) divergence, between the true distributions $A_j \lambda$ and their empirical estimates $p_j$.


By \eqref{eqn:pq_{ij}}, each item-assortment pair $i \in \CA_j$ is seen $Kq_{ij}$ times amongst the observations $\left\{i^k,\cA^k\right\}_{k=1}^K$. Based on this, the log-likelihood of the observation set $\left\{i^k,\cA^k\right\}_{k=1}^K$ is
$\sum_{j \in [m]} \sum_{i \in \cA_j} Kq_{ij} \log\left( \la a_{ij}, \lambda \ra\right)$. 
Thus, ignoring the constant $K$ factor, the MLE problem is
\begin{align}\label{eqn:MLE-frequency}
\max_\lambda \left\{ \sum_{j \in [m]} \sum_{i \in \cA_j} q_{ij} \log\left( \la a_{ij}, \lambda \ra\right) :~  \lambda \in \Delta_{n!} \right\}.
\end{align}
Throughout, we use the convention that when $q_{ij} = \la a_{ij}, \lambda \ra= 0$, we set $q_{ij} \log(\la a_{ij}, \lambda\ra) = 0$. This implies that if the optimal solution $\lambda$ to \eqref{eqn:MLE-frequency} has $\bbP_\lambda[i \mid \cA_j] = \la a_{ij}, \lambda \ra= 0$, then we must have $q_{ij} = 0$ also, i.e., we did not observe any choices of $i$ from $\cA_j$ in our data either.

Like \eqref{eqn:revenue-prediction-model}, the problem \eqref{eqn:MLE-frequency} is very large, with $n!$ variables. A column generation technique is suggested
in \citet{vanRyzinVulcano2015} 
to get around this, i.e., solve \eqref{eqn:MLE-frequency} on a subset of the variables, and use the optimality conditions to add variables as needed.
The MLE column generating subproblem is constructed as
\begin{equation}\label{eqn:MLE-column-generating}
\max_{\sigma} \left\{ \sum_{j \in [m]} \sum_{i \in \cA_j} \frac{q_{ij} a_{ij}(\sigma)}{\la a_{ij}, \lambda(S)\ra} :~ \sigma \in S_n \right\}.
\end{equation}
The solution $\lambda(S)$ is optimal if $\eqref{eqn:MLE-column-generating} \leq K$, otherwise the column $\sigma^*$ maximizing \eqref{eqn:MLE-column-generating} is added to the set $S$, and the process is repeated. 
Note that \eqref{eqn:MLE-column-generating} has the same form as \eqref{eqn:subgradient-subproblem} and \eqref{eqn:revenue-prediction-robust-constraint}.

We next demonstrate that the MLE problem \eqref{eqn:MLE-frequency} admits a nice interpretation between the empirical estimates $\left\{ p_j \right\}_{j \in [m]}$ and the distributions $\left\{ A_j \lambda \right\}_{j \in [m]}$. To observe this, let us rewrite the objective in \eqref{eqn:MLE-frequency} as
\begin{align*}
\sum_{j \in [m]} \sum_{i \in \cA_j} q_{ij} \log\left(\la a_{ij}, \lambda \ra \right) 
&=  \sum_{j \in [m]} q_j \sum_{i \in \cA_j} p_{ij} \log\left( \la a_{ij}, \lambda \ra  \right)\\
&=  - \sum_{j \in [m]} q_j \underbrace{\sum_{i \in \cA_j} p_{ij} \log\left( \frac{p_{ij}}{\la a_{ij}, \lambda \ra } \right)}_{=  \KL(p_j,A_j \lambda)}  + \underbrace{\sum_{j \in [m]} q_j \sum_{i \in \cA_j} p_{ij} \log(p_{ij})}_{= \text{constant}}
\end{align*}
where $\KL(a,b)$ is the KL divergence between two probability distributions $a$ and $b$.
Hence, \eqref{eqn:MLE-frequency} is equivalent to solving
\begin{equation}\label{eqn:MLE-frequency-alt}
\min_{\lambda} \left\{ \sum_{j \in [m]} q_j \KL(p_j,A_j \lambda) :~  \lambda \in \Delta_{n!} \right\}.
\end{equation}
Thus, by defining $D(A\lambda,p) = \sum_{j \in [m]} q_j \KL(p_j,A_j \lambda)$, we see that the MLE approach is equivalent to \eqref{eqn:revenue-prediction-distance-min} but with a different distance metric $D(\cdot,\cdot)$.

\subsection{Norm-Minimization Approach}\label{sec:norm-min-approach}

As opposed to the approaches outlined in Sections~\ref{sec:rev-predict}~and~\ref{sec:MLE-approach}, in order to estimate a non-parametric choice model $\lambda$, \citet{BertsimasMisic2015} suggest minimizing the $\ell_1$-norm of $p - A \lambda$ by solving 
\begin{equation}\label{eqn:norm-l1-min}
\min_{\lambda} \left\{ \|p - A \lambda\|_1 :~  \lambda \in \Delta_{n!} \right\}.
\end{equation}


In fact, \eqref{eqn:norm-l1-min} can be cast as an LP, but it is still computationally intractable since the dimension of $\lambda$ is $n!$. Similar to \citet{vanRyzinVulcano2015}, \citet{BertsimasMisic2015} addresses this computational difficulty via a column generation approach. Again, \eqref{eqn:norm-l1-min} is of the same form as \eqref{eqn:revenue-prediction-distance-min} where the distance metric $D(\cdot,\cdot)$ is selected to be $D(A\lambda,p) = \|p - A\lambda\|_1$.
Furthermore, the resulting column generating subproblem is of the form
\begin{equation}\label{eqn:norm-l1-min-column-generating}
\max_{\sigma} \left\{ \sum_{j \in [m]} \sum_{i \in \cA_j} \beta_{ij}(S) a_{ij}(\sigma) - \nu(S) : ~ \sigma \in S_n \right\},
\end{equation}
where $\beta(S)$ and $\nu(S)$ are from the dual solution to solving \eqref{eqn:norm-l1-min} on a subset of columns $\sigma \in S \subset S_n$. Again, this subproblem has the same form as \eqref{eqn:subgradient-subproblem}, \eqref{eqn:revenue-prediction-robust-constraint} and \eqref{eqn:MLE-column-generating}.





Experiments.tex
\section{Computational Study}\label{sec:experiments}

In our computational study, we compare three different choices of distance functions $D(\cdot,\cdot)$ arising from three different $\ell_p$-norms $\ell_1,\ell_2,\ell_\infty$, as well as some parametric MNL approaches. In the static data setting, we examine model fit for all approaches, and sparsity properties for the non-parametric approaches. We also examine model fit and sparsity properties of the non-parametric approach in the dynamic data setting. Because the convergence properties of the dual approach (both in terms of different choices of $D$ and also the dynamic setup) are superior than those of the primal approach, in these experiments we only employ the dual approach. All experiments are conducted on a server with 2.8 GHz processor and 64GB memory, using Python 3.5. Gurobi 7.0 is used to solve the integer programming subproblems.

We employ the same simulated data setup from \citet[Section 5.3]{BertsimasMisic2015}. Our ground truth choice model is a mixed MNL model generated by two meta-parameters: $K \in \bbN$, the number of MNL models in the mixture, and $L \in \bbR_+$ governing the intensity of the preferences. Given mixing probabilities $w \in \Delta_K$ and $K$ sets of utilities $\{u_{i,k}\}_{i \in \{0\} \cup [n]}$, $k \in [K]$, the mixed MNL model chooses an item $i \in \CA \subseteq [n]$ with probability
\vspace{-7pt}
\[ \bbP[i \mid \CA] = \sum_{k \in [K]} w_k \frac{u_{i,k}}{u_{0,k} + \sum_{i' \in \CA} u_{i',k}}.\]
\vspace{-13pt}

For each $k \in [K]$, we generate $n+1$ parameters $q_{i,k} \sim U(0,1)$, $i \in \{0\} \cup [n]$ (recall that $0$ denotes the no-choice option present in each subset). The utilities $u_{i,k}$ are then set as follows: four randomly chosen $i \in \{0\} \cup [n]$ are set to $u_{i,k} = Lq_{i,k}$, while the rest are set to $u_{i,k} = q_{i,k}/10$. The mixing probabilities $\{w_k\}_{k \in [K]}$ are chosen randomly from the $(K-1)$-dimensional simplex. We test the case $n=10$ and parameter regimes $K \in \{1,5,10\}$ and $L \in \{5,10,100\}$.

For each combination of $K$ and $L$, we generate 100 instances of the ground truth mixed MNL model. For each instance we generate 150 subsets of maximum size $\lfloor n/2 \rfloor$ uniformly at random. We reserve 50 subsets for training, and the other 100 for testing. Our training regime is as follows: using the ground truth model, we compute the $p_{\train}$ vector, where $p_{\train,ij} = \bbP[i \mid \CA_j]$, and $\CA_j$ is are subsets from our training set. We vary the number of training subsets $m$ we use by $m \in \{10,20,50\}$.

We use $p_{\train}$ to fit our model using the methods from Section \ref{sec:learning-choice-model-MD}. To evaluate model fit, we similarly compute $p_{\test}$ for our ground truth model using all 100 subsets in the test set, then compute $\hat{p}_{\test}$ from our learned model, and examine the mean average error (MAE) between $p_{\test}$ and $\hat{p}_{\test}$, defined as $\MAE(p, \hat{p}) = \frac{1}{\text{length}(p)} \sum_{i,j} |p_{ij} - \hat{p}_{ij}|$.
We terminate training when $\MAE(p_{\train}, \hat{p}_{\train}) \leq 0.001$, and we evaluate performance by computing $\MAE(p_{\test}, \hat{p}_{\test})$.

\input{comparefig1}

Figure \ref{fig:MAE-m_vs_p} displays the test set MAE for different $K$, $m$ and $\ell_p$ while $L=5$ is fixed. As expected, we see that as $m$ increases, the test MAE decreases. This trend holds for all combinations of $L$, $K$ and $\ell_p$, so we next examine results for for $m=20$.
\input{comparefig2}

Figure \ref{fig:MAE-m_vs_p} suggests a slight upward trend as we increase our $\ell_p$-norm. To study this further, we plot the average test MAE and the average number of rankings in the estimated choice model at termination across the 100 instances for $m=20$ and each parameter combination of $K$ and $L$ in Figure \ref{fig:MAE_numrankings_compare}. While the scales of the left $y$-axes are quite small, and hence the differences are quite minor, we notice that $\ell_1$ often has the best test MAE. The two exceptions are Figures \ref{fig:compare_L10K1} and \ref{fig:compare_L100K1}, where $\ell_2$ has the best fit. Also notice that in all cases except Figure \ref{fig:compare_L100K1}, $\ell_\infty$ generally has worse fit than $\ell_1,\ell_2$.
While $\ell_1$ often fits better than $\ell_2,\ell_\infty$, it comes at the cost of a significantly denser final model. Figure \ref{fig:MAE_numrankings_compare} shows that often the number of rankings for $\ell_1$ is more than twice that of $\ell_2,\ell_\infty$, which are generally quite similar.
Thus, choosing $\ell_1$ will give us the best fitting model on our test set, but comes at the cost of a significant loss of sparsity. Instead, $\ell_2$ seems to lead to a good compromise, for which we observe similar test MAE as $\ell_1$, but generally learns a sparser model.

Next, we compare our non-parametric choice model estimation technique to fitting a parametric mixed MNL model trained using maximum likelihood estimation \citep{Train2009}. This involves a selection of hyperparameter $K_{\fit}$, the number of MNL segments. We test $K_{\fit} \in \{1,3,5,10\}$. Figure \ref{fig:MAE-par_vs_nonpar} compares the test MAE of the estimated models for different $m$ and $K$, holding $L=5$ fixed.  

\input{comparefig3}

As we expect, choosing higher $K_{\fit}$ reduces the MAE for parametric models, and higher ground truth $K$ has a more pronounced effect on the parametric models than the non-parametric ones. However, notice that the most dramatic effect is by changing $m$, the number of training sets. For $m=10$, the parametric models outperform the non-parametric models; for $m=20$, the two types of models are roughly even; but for $m=50$, the non-parametric model outperforms the parametric ones. This suggests that non-parametric models are better at exploiting settings where we have an abundance of data.

We also examine our non-parametric approach in the dynamic data setting when we only have access to $p^t$. We fix the norm to be $\ell_2$-norm, $L=5$, $K=5$, and $m=20$. We generate $p^1$ by randomly generating $2000$ observations of item-assortment choices from the ground truth model, and aggregating them according to \eqref{eqn:pq_{ij}}. Each subsequent $p^t$ is generated by adding $\kappa \in \bbN$ observations to the previous $p^{t-1}$. We stop when $\MAE(\bar{p}^T, \hat{p}) \leq 0.001$, where $\bar{p}^T$ is the average of the $p^t$ vectors seen. We test the effects of using different $\kappa \in \{50,100,200,300,400,500,600,700,800,900,1000,\infty\}$, where $\kappa = \infty$ denotes using the true vector $p$.
\input{comparefig4}

Figure \ref{fig:dynamic-run} shows that the average test MAE, once the stopping criterion is reached, are virtually identical for different $\kappa$. Figure \ref{fig:dynamic-increments} plots the number of iterations for different $\kappa$, and shows that the learned model becomes sparser when we increase $\kappa$. Note that the runs for $\kappa=50$ take roughly twice as many iterations as $\kappa=1000$, which means that runs for $\kappa = 1000$ use roughly $10$ times as many observations to achieve effectively the same test MAE error. Surprisingly, using the \emph{true} probabilities $p$ to learn a model also achieves the same test MAE. We deduce that the main gains to additional observations are faster convergence and sparsity of our learned model, but test set accuracy remains unaffected. Also, the gains to sparsity diminish rapidly as $\kappa$ increases.





Intro.tex
\section{Introduction}\label{sec:intro}

A choice model specifies the probability distribution for rankings that consumers have over a set of items. Such models give choice probabilities, that is, the probability that a given consumer will choose an item from a given subset. Choice models are prevalent in several application areas such as revenue management, web page ranking, betting theory, social choice, marketing, and economics (see \citet{DworkKumarNaorSivakumar2001,TallurivanRyzin2006book,JagabathulaShah2008,FariasJagShah2009,DesirGoyalJagSegev2016} and references therein). A good choice model aims to capture complex substitution behaviors of consumers in order to accurately predict demand from limited observations.

Choice model estimation has received quite a bit of interest. Traditional choice models often specify a \emph{parametric} structure for the probability distribution (examples include the multinomial logit (MNL),
nested logit, and mixed MNL models), see \citet{TallurivanRyzin2006book}
and references therein.
Imposing a parametric structure makes estimation of the necessary parameters a simpler task, but is often at the expense of 
overly facile assumptions on consumer behavior (such as independence of irrelevant alternatives in MNL models) preventing us from
accurately capturing the substitution behavior.
Because of this, the \emph{non-parametric} approach of directly estimating the probability distribution for rankings have drawn growing interest in academia and practice \citep{Rusmevichientong2006,FariasJS2013},  
and is shown through case studies \citep{FariasJS2013,JagabathulaRusmevichientong2016}
to lead to substantial improvement in prediction accuracy.
However, learning (or even specifying) a full non-parametric model is intractable even for moderate-sized problems, since there are a factorial number of probabilities to estimate.
Furthermore, today's technology enables us to collect data on a continuous basis. Unfortunately, existing methods have a notable deficiency: in order to update their estimates with the latest data, they must \emph{re-solve} their choice models \emph{entirely}. In the case of non-parametric choice estimation, such dynamic considerations significantly compound the existing computational challenges.


\textbf{Related literature.} Earliest studies on non-parametric choice models appear in the economics and psychology literatures, e.g., \citet{BlockMarschak1960}. \citet{MahajanVanRyzin2001} showed that non-parametric models capture a number of parametric models as special cases.
There are three main approaches 
employed for the static estimation of non-parametric choice models. These approaches are essentially characterized with different objective functions.  
\citet{FariasJS2013} outline a constraint sampling based method to estimate the non-parametric model from observational data on consumer choices. Other methods based on maximum likelihood estimation (MLE) and norm minimization are presented in \citet{vanRyzinVulcano2015} and \citet{BertsimasMisic2015} respectively.
\citet{JagabathulaRusmevichientong2016} suggest to use the Frank-Wolfe algorithm to estimate a non-parametric choice model with a general distance function, which captures previous approaches \citep{vanRyzinVulcano2015,BertsimasMisic2015}.  
While \citet{FariasJS2013} and \citet{vanRyzinVulcano2015} provide useful recovery results under some assumptions on the observational data, none of the methods in \citet{FariasJS2013, vanRyzinVulcano2015, BertsimasMisic2015, JagabathulaRusmevichientong2016} come with provable convergence guarantees. Moreover, these prior methods are not equipped to deal with the more realistic \emph{dynamic} data setting, where the firm continuously collect and wishes to utilize more consumer data even as we implement an estimation procedure.




\textbf{Our contributions.} We present two simple iterative frameworks for estimating non-parametric choice models via minimizing a general distance function in the usual \emph{static} setting as well as the \emph{dynamic} setting, where new empirical observations are continuously collected. We complement the existing literature by providing explicit error bounds (in the form of convergence rates) for the recovery performance of our iterative approaches in both static and dynamic settings.

Our first approach follows a primal approach to the minimization problem via the Frank-Wolfe (F-W) algorithm. This approach was recently suggested in \citet{JagabathulaRusmevichientong2016} for the static estimation problem. We show that, even in the static setup, for two classes of distance functions suggested previously in \citet{BertsimasMisic2015,vanRyzinVulcano2015}, the F-W algorithm does not enjoy the usual convergence rate described in \citet{Jaggi2013} due to the unboundedness of an important quantity (the curvature constant) associated with the standard analysis of the F-W algorithm. 
Nonetheless, we identify a set of assumptions that ensure efficient convergence guarantees for this method. In particular, we establish that our convergence guarantees for the F-W algorithm apply to a class of distance functions based on norms. Furthermore, we show that the F-W method adapts naturally to the \emph{dynamic} setting. In this latter setting, we establish an error bound on a dynamic variant of F-W algorithm that depends on the rate at which the updated empirical choice probabilities converge to a limit. We demonstrate that convergence of this method in the dynamic setup is achievable as long as the data converges at a certain minimum rate.

To alleviate the need for a minimum rate of convergence of data in the case of F-W algorithm, we analyze a  different approach. Our second approach examines the non-parametric estimation problem via duality, convex conjugacy and regret minimization from online learning. Again, we prove convergence of this method in the dynamic setting by providing efficient guarantees on the number of iterations, but with no dependence or requirement on the minimum rate at which data converges. 

Both our primal and dual approaches achieve the same convergence rate in the static setup in spite of different assumptions used in their dynamic versions. Moreover, in both of our approaches, by establishing the number of iterations needed in order to estimate the choice model to within a given accuracy $\epsilon$, our analysis also upper-bounds the sparsity of our estimated model, and as a result, it exposes an explicit trade-off between the simplicity of the non-parametric choice model and its estimation accuracy.


We also carry out a computational study, where 
we test the behavior of different distance metrics ($\ell_p$-norms) in terms of estimation accuracy and sparsity of our learned models as well as the impact of number of observations and batch size used in the dynamic data setting on the same criteria. 
Our results indicate that several choices in designing a generic dynamic setup mainly influence the simplicity of the estimated choice model as opposed to its accuracy.



\textbf{Notation.} For a positive integer $n\in\N$,  we let $[n]=\{1,\ldots,n\}$, define $\Delta_n:=\{x\in\R^n_+:~\sum_{i\in[n]} x_i=1\}$ to be the standard simplex, and $S_n$ to be the collection of rankings of the set $[n]$.
We refer to a collection of objects $b_j$, $j \in J$ by the notation $\{b_j\}_{j \in J}$. Throughout the paper, the superscript, e.g., $y^t,z^t,f^t$, is used to attribute items to the $t$-th time period or iteration. 
The subscript is used to denote coordinates of a vector or matrix, e.g., $\beta_{ij}$. Given vectors $x$ and $y$, $\la x,y\ra$ corresponds to the usual inner product of $x$ and $y$. Given a norm $\|\cdot\|$ on a Euclidean space $\bbE$ and a real number $a> 0$, we
denote its dual norm by $\|x\|_* = \min_{y} \{ \la x, y\ra:~\|x\| \leq1 \}$.
For $q\in[1,\infty]$, $\|x\|_q$ denotes the $\ell_q$ norm of $x$. 
We let $\partial f(x)$ be the subdifferential of $f$ taken at $x$.
We abuse notation slightly by denoting $\grad f(x)$ for both the gradient of function $f$ at $x$ if $f$ is differentiable and a subgradient of $f$ at $x$, even if $f$ is not differentiable. If $\phi$ is of the form $\phi(x,y)$, then $\grad_x \phi(x,y)$ denotes the subgradient of $\phi$ at $x$ while keeping the other variables fixed at $y$.
We denote the indicator function as $\I$, i.e., $\I(\cS) = 1$ if statement $\cS$ holds, and $\I(\cS) = 0$ otherwise.


\section{Data and Model}\label{sec:model-data}



In the general non-parametric choice estimation framework, we consider a firm which has $n$ products to sell. We assume that the firm chooses from a given set of $m$ assortments $\cA_1,\ldots,\cA_m \subset [n]$, and the item $1$ represents the `no-buy' option and is present in all assortments $\cA_j$, $j \in [m]$. Therefore, when presented with an assortment, the consumer will always choose an item from it (perhaps the no-buy option, i.e., item $1$). We assume for simplicity that the order of items in an assortment, their duration of appearance and the number of times an item appears in different assortments play no role in the choice. We also denote $N := \sum_{j=1}^{m} |\cA_j|$.

In this paper, we work with the following data collection process. When a consumer arrives, the firm displays an assortment $\cA_j$ to the consumer. The consumer chooses a product $i \in \cA_j$, and the firm observes this choice.
Our static data set consists of $K$ such observations, i.e., a collection of pairs $\left\{i^k,\cA^k\right\}_{k=1}^K$, where $i^k$ denotes the item chosen and $\cA^k$ denotes the assortment displayed for observation $k$. There are a number of useful statistics on this data set, which are defined as follows:
\begin{equation}
q_{ij} := \frac{1}{K} \sum_{k=1}^{K} \I(i^k = i, \cA^k = \cA_j), \quad q_{j} := \frac{1}{K} \sum_{k=1}^{K} \I(\cA^k = \cA_j) \quad\text{and}\quad p_{ij} := \frac{q_{ij}}{q_j}.\label{eqn:pq_{ij}}
\end{equation}
In words, $q_{ij}$ is the proportion of observations where assortment $\cA_j$ was displayed and item $i$ was chosen, $q_j$ is the proportion of observations where assortment $\cA_j$ was displayed, and $p_{ij}$ is the proportion of consumers who chose item $i$ \emph{given that} assortment $\cA_j$ was displayed.

The non-parametric choice model is as follows. An incoming consumer will choose an item according to his/her ranking $\sigma \in S_n$ of the products in $[n]$, that is, when presented with an assortment $\cA_j$, the consumer chooses the highest ranked product from $\cA_j$, i.e., $\argmin_{i \in \cA_j} \sigma(i)$. The key assumption in this model is that the ranking $\sigma$ of each incoming consumer is distributed i.i.d.\@ according to some distribution $\lambda$ on the set of all rankings $S_n$. Here $\lambda$ represents the vector of probabilities of each $\sigma \in S_n$ being drawn, i.e., the probability that a consumer will have a ranking $\sigma$ is $\lambda(\sigma)$.





For an item-assortment pair $i \in \cA_j$ and ranking $\sigma \in S_n$, we define $a_{ij}(\sigma) = 1$ if $i$ is the highest ranked item in $\cA_j$ according to $\sigma$, and $a_{ij}(\sigma) = 0$ otherwise. We define $A$ to be the binary matrix of dimension $N \times n!$ with entries $a_{ij}(\sigma)$. Each row corresponds to a pair $i \in \CA_j$ and each column corresponds to a ranking $\sigma$. We denote the rows as $a_{ij}^\top \in \{0,1\}^{n!}$ and the columns as $a(\sigma) \in \{0,1\}^N$.
For $j \in [m]$, we define $A_j$ to be the submatrix of $A$ with rows $a_{ij}^\top$ for $i \in \cA_j$. 

Based on our notation, given a distribution $\lambda \in \Delta_{n!}$ over rankings, the choice probability $\bbP_{\lambda}[i \mid \cA_j]$ that a random consumer chooses item $i$ from $\cA_j$ is simply represented as the inner product $\la a_{ij}, \lambda\ra$. Moreover, $A \lambda$ is simply the vector composed of the choice probabilities $\left\{\bbP_{\lambda}[i \mid \cA_j]\right\}_{i \in \cA_j, j \in [m]}$. The choice probabilities $\left\{ \bbP_\lambda[i \mid \cA_j] \right\}_{i \in \cA_j}$ for each $\CA_j$ is simply $A_j \lambda$.
Finally, we define the polytope of all possible choice probabilities consistent with a distribution as
\[X := \left\{ A\lambda :~ \lambda \in \Delta_{n!} \right\} = \conv\left( \left\{a(\sigma) :~ \sigma \in S_n \right\} \right) \subseteq \bbR^N.\]
The statistics \eqref{eqn:pq_{ij}} will be used to infer the best-fitting probability distribution $\lambda$, since $p_{ij}$ is an empirical estimate of $\bbP_\lambda[i \mid \CA_j]$. We denote the collected vectors of $\left\{p_{ij}\right\}_{i \in \cA_j, j \in [m]} = p \in \bbR^N$ and $\left\{p_{ij}\right\}_{i \in \cA_j} = p_j \in \bbR^{|\cA_j|}$ for $j \in [m]$.






NonparChoice.tex
\documentclass[11pt,oneside]{article}

\usepackage{amsthm,amsmath,amssymb,amsfonts}
\usepackage{mathtools,xfrac}
\usepackage{url,setspace,color}
\usepackage{graphicx,epsfig}
\usepackage{enumitem}
\usepackage{pdfsync,array}
\usepackage[numbers,sort&compress]{natbib}
\usepackage[normalem]{ulem}
\usepackage{algorithm,algorithmic}
\newcommand{\LINEIF}[2]{%
    \STATE\algorithmicif\ {#1}\ \algorithmicthen\ {#2}%
}

\usepackage{subfigure,caption}
\usepackage{bm}
\usepackage{tikz}
\usepackage{pgfplots}
\pgfplotsset{compat=1.3}

\DeclareMathOperator{\train}{train}
\DeclareMathOperator{\test}{test}
\DeclareMathOperator{\MAE}{MAE}


\usepackage{geometry,verbatim}
\usepackage{makeidx,latexsym}
\usepackage{bm}
\usepackage{booktabs}
\usepackage{multirow}
\usepackage{subfigure}

\usepackage[colorlinks=true]{hyperref}

\usepackage{authblk}
\usepackage[in]{fullpage}
\makeatletter
\def\imod#1{\allowbreak\mkern10mu({\operator@font mod}\,\,#1)}
\makeatother
\usepackage[compact]{titlesec}

\makeatletter
\renewenvironment{proof}[2][\proofname]{\par
\pushQED{\qed}
\normalfont
\topsep6\p@\@plus6\p@\relax
\trivlist
\item[\hskip\labelsep\textit{#2}\@addpunct{.}]\ignorespaces
}{%
\popQED\endtrivlist\@endpefalse
}
\makeatother


\def\halmos{\mbox{\quad$\square$}}
\def\Halmos{\mbox{\quad$\square$}}
\def\Halmos{\qed}
\def\Halmos{}


\theoremstyle{plain}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{conjecture}{Conjecture}[section]
\newtheorem{lemma}{Lemma}[section]
\newtheorem*{lemma*}{Lemma}
\newtheorem{proposition}{Proposition}[section]
\newtheorem*{proposition*}{Proposition}
\newtheorem{corollary}{Corollary}[section]
\newtheorem{observation}{Observation}[section]
\newtheorem*{observation*}{Observation}
\newtheorem{assumption}{Assumption}[section]
\newtheorem{condition}{Condition}[section]


\theoremstyle{definition}
\newtheorem{question}{Question}[section]
\newtheorem{definition}{Definition}[section]

\theoremstyle{remark}
\newtheorem{remark}{Remark}[section]
\newtheorem*{remark*}{Remark}
\newtheorem{example}{Example}[section]
\newtheorem*{example*}{Example}




\def\VRHDW#1#2#3{\vrule height #1 depth #2 width #3}
\newcommand{\up}{\VRHDW{1.5em}{0em}{0em}}
\newcommand{\down}{\VRHDW{0em}{0.625em}{0em}}%
\def\updown{\up\down}
\def\TABLE{\caption}


\input{Def}
\usepackage[colorlinks=true]{hyperref}
\def\crd{\color{black}}
\def\cbl{\color{black}}
\def\cgr{\color{black}}
\def\cmg{\color{black}}
\def\ccy{\color{black}}
\def\cye{\color{black}}


\begin{document}

\title{Dynamic Data-Driven Estimation of Non-Parametric Choice Models}
\author[1]{Nam Ho-Nguyen}
\author[1]{Fatma K{\i}l{\i}n\c{c}-Karzan}
\affil[1]{Tepper School of Business, Carnegie Mellon University, Pittsburgh, PA, 15213, USA.}
\date{February 13, 2017; latest version February 21, 2018}
\maketitle

\begin{abstract}
We study non-parametric estimation of consumer choice models. Non-parametric approaches were introduced to alleviate unreasonable assumptions and issues of suboptimal model fit/selection present in traditional parametric approaches, and are prevalent in several application areas. We present two convex optimization-based frameworks, a primal approach and a dual approach, to efficiently learn a non-parametric choice model from data that is close to the best-fitting one. As opposed to the existing literature, both approaches enjoy provable convergence guarantees and extend naturally to the \emph{dynamic} observation setting. Our computational study on the dynamic setting reveals the true impact of how much data are needed and at what rate to achieve the best trade-off in terms of estimation accuracy and model simplicity. In the static setup, we also compare our non-parametric approach with existing parametric approaches.
\end{abstract}

\maketitle

\input{Intro}

\input{DynamicLearning}

\input{Experiments}


\section*{Acknowledgements}

This research is supported in part by NSF grant CMMI 1454548.

\bibliographystyle{abbrvnat}
\input{NonparChoice.bbl}

\newpage
\appendix

\input{ExistingApproaches}
\input{Proofs}
\input{OCOframework}



\end{document}






OCOframework.tex
\section{Online Convex Optimization Framework}\label{sec:prox-setup}

Online convex optimization (OCO) is commonly used to capture decision making in dynamic environments. Here we outline the basic OCO concepts; for further details and background, we refer the reader to \citet{CesaBianchiLugosi2006book}. 

In OCO, we are given a finite time horizon $T$, closed, bounded, and convex domain $Z$, and in each time period $t \in [T]$, a convex loss function $f^t:Z \to \R$ is revealed. At time periods $t \in [T]$ we must choose a decision $z^t \in Z$, and based on this we suffer a loss of $f^t(z^t)$ and receive some feedback typically in the form of first-order information on $f^t$. The main aim of OCO is to choose a sequence of points $\{z^t\}_{t=1}^T$ from the domain $\cZ$ to bound the weighted regret
\begin{equation}\label{eqn:weighted-regret}
\sum_{t=1}^{T} \theta^t f^t(z^t) - \inf_{z \in \cZ} \sum_{t=1}^{T} \theta^t f^t(z),
\end{equation}
where $\theta = \{\theta^t\}_{t=1}^T \in \Delta_T$ is a collection of convex combination weights. The key restriction that separates OCO from standard optimization problems is that $z^t$ must be chosen \emph{before} observing $f^t$. The fact that there exist algorithms which bound \eqref{eqn:weighted-regret} for \emph{any} sequence $\{f^t\}_{t=1}^T$ is the crucial  aspect of OCO which we exploit to solve the dynamic variant of \eqref{eqn:norm-min-convex-form}.

A key class of algorithms which can be used for OCO (as well as standard offline convex optimization) are first-order methods (FOMs). Following the notation in the survey of \citet{JuditNem2012Pt1}, we outline the proximal setup for a general domain $\cZ$. This setup forms the basis for several FOMs such as Mirror Descent and is used in their convergence analyses.
\begin{itemize}
\item \emph{Norm}: $\|\cdot\|$ on the Euclidean space $\E$ where the domain $\cZ$ lives, along with its dual norm $\|\zeta\|_*:=\max\limits_{\|z\|\leq1}\langle\zeta,z\rangle$.
\item \emph{Distance-Generating Function} ({\dgf}): A function $\omega(z):\cZ\rightarrow \R$, which is convex and continuous on $\cZ$, and admits a selection of subgradients  $\grad\omega(z)$ that is continuous on the set $\cZ^\circ:=\{z\in \cZ:\partial\omega(z)\neq\emptyset\}$ (here $\partial\omega(z)$ is a subdifferential of $\omega$ taken at $z$), and is strongly convex with modulus 1 with respect to  
$\|\cdot\|$:
\[
\forall z',z''\in \cZ^\circ:~ \langle \grad\omega(z')-\grad\omega(z''),\ z'-z''\rangle \geq\|z'-z''\|^2. 
\]
\item  \emph{Prox-mapping}: Given a \emph{prox center} $z\in \cZ^\circ$,
\[
\Prox_z(\xi):=\argmin\limits_{z' \in \cZ}\left\{\langle \xi-\grad\omega(z),z'\rangle + \omega(z') \right\}: \E\to \cZ^\circ.
\]
When the {\dgf} is taken as the squared $\ell_2$-norm, the prox mapping becomes the usual projection operation of the vector $z-\xi$ onto $\cZ$.
\item \emph{Set width}: $\Omega=\Omega_z := \max\limits_{z\in \cZ}\omega(z)-\min\limits_{z\in \cZ}\omega(z)$.
\end{itemize}
For common domains $\cZ$ such as simplex, Euclidean ball, and spectahedron, standard proximal setups, i.e., selection of norm $\|\cdot\|$, {\dgf} $\omega(\cdot)$, the resulting $\Prox$ computations and set widths $\Omega$ are discussed in \cite[Section 1.7]{JuditNem2012Pt1}. 

In the most basic setup, our functions $f^t$ are convex and non-smooth. In this case, we utilize a generalization of Mirror Descent, outlined in Algorithm~\ref{alg:md} for bounding the weighted regret  \eqref{eqn:weighted-regret}.
\begin{algorithm}[t]
\caption{Generalized Mirror Descent}\label{alg:md}
\begin{algorithmic}
\STATE {\bf Input:~} time horizon $T$, positive step sizes $\{\gamma^t\}_{t=1}^T$, and a sequence of vectors $\{\xi^t\}_{t=1}^T$.
\STATE {\bf Output:~} sequence $\{z^t\}_{t=1}^T$ from $\cZ$.
\STATE $z^1 := \min_{z \in \cZ} \omega(z)$.
\FOR{$t=1,\ldots,T$}
	\STATE $z^{t+1} = \Prox_{z^t} \left( \gamma^t \xi^t \right)$.
\ENDFOR
\end{algorithmic}
\end{algorithm}

We next state a bound on the weighted regret \eqref{eqn:weighted-regret} in the most general case where our functions $f^t$ need only satisfy convexity and Lipschitz continuity. More precisely, we will assume the following.
\begin{assumption}\label{ass:OCO-non-smooth}\label{cond:prox-setup}
A proximal setup of Section~\ref{sec:prox-setup} exists for the domain $\cZ$. Each function $f^t$ is convex, and there exists $G \in(0,\infty)$ such that the subgradients of  $f^t$ are bounded, i.e., $\|\grad f^t(z)\|_* \leq G$ for all $z \in \cZ$ and $t \in [T]$.
\end{assumption}

\begin{theorem}[{\citet[Theorem 1]{HoNguyenKK2016}}]\label{thm:OCO-non-smooth}
Suppose Assumption~\ref{ass:OCO-non-smooth} holds, and we are given weights $\theta \in \Delta_T$. Then
Algorithm~\ref{alg:md} with $\xi^t = \theta^t \grad f^t (z^t)$, and step sizes $\gamma_t =\gamma:= \sqrt{\frac{2 \Omega}{\sup_{t \in [T]} (\theta^t)^2 G^2 T}}$ for all $t\in[T]$ results in
\[ \sum_{t=1}^{T} \theta^t f^t(z^t) - \inf_{z \in \cZ} \sum_{t=1}^{T} \theta^t f^t(z) \leq \sqrt{2 \Omega \left(\sup_{t \in [T]} (\theta^t)^2\right) G^2 T}. \]
\end{theorem}
The bound on weighted regret in Theorem~\ref{thm:OCO-non-smooth} is optimized when the convex combination weights $\theta \in \Delta_T$ are set to be \emph{uniform}, i.e., $\theta^t = 1/T$; in this case, the bound
above becomes $O(1/\sqrt{T})$.

Proofs.tex
\section{Proofs for Section \ref{sec:learning-choice-model-MD}}\label{sec:proofs}

\proof{Proof of Theorem~\ref{thm:FW-dynamic-convergence}.}



From Assumption \ref{ass:distance-smoothness} applied to $D(\cdot,p^t)$ and $x^{t+1} = (1-\gamma^t)x^t + \gamma^t z^t$, we get
\[D(x^{t+1},p^t) \leq D(x^t,p^t) + \gamma_t \la \grad_x D(x^t,p^t), z^t - x^t \ra + (\gamma^t)^2 C_D. \]
From the definition of $z^t$ and convexity of $D(\cdot,p^t)$, we have $\la \grad_x D(x^t,p^t), z^t - x^t \ra \leq \la \grad_x D(x^t,p^t), x - x^t \ra \leq D(x,p^t) - D(x^t,p^t)$ for any $x \in X$, so
\[D(x^{t+1},p^t) \leq (1-\gamma_t) D(x^t,p^t) + \gamma_t D(x,p^t) + (\gamma^t)^2 C_D. \]
Adding and subtracting the appropriate terms, we arrive at
\begin{align*}
&D(x^{t+1},p) - D(x,p)\\
&\leq (1-\gamma_t) (D(x^t,p) - D(x,p)) + (\gamma^t)^2 C_D\\
&\quad + D(x^{t+1},p) - D(x^{t+1},p^t) + (1-\gamma_t)(D(x^t,p^t) - D(x^t,p))\\
&\quad + \gamma_t (D(x,p^t)-D(x,p)).
\end{align*}
Defining $\delta_t := D(x^t,p) - D(x,p)$ and $\alpha_t := (\gamma^t)^2 C_D + D(x^{t+1},p) - D(x^{t+1},p^t) + (1-\gamma_t)(D(x^t,p^t) - D(x^t,p)) + \gamma_t (D(x,p^t)-D(x,p))$, we now have the recursion $\delta_{t+1} \leq (1-\gamma_t) \delta_t + \alpha_t$.

By induction, for $t \geq 3$,
\[\delta_t \leq \delta_1 \prod_{k=1}^{t-1}(1-\gamma_k) + \sum_{k=1}^{t-2} \alpha_k \prod_{l=k+1}^{t-1} (1-\gamma_l) + \alpha_{t-1}.\]
Moreover, the first term disappears because $\gamma_k = 2/(k+1)$, $1-\gamma_1 = 0$. Furthermore, $\prod_{l=k+1}^{t-1} (1-\gamma_l) = k(k+1)/((t-1)t)$ for $k=1,\ldots,t-3$ and $(t-2)/t$ for $k=t-2$, so
\[\delta_t \leq \frac{1}{(t-1)t} \sum_{k=1}^{t-3} k(k+1) \alpha_k + \frac{t-2}{t} \alpha_{t-2} + \alpha_{t-1}.\]
Substituting the definitions of $\delta_t,\alpha_t$ and simplifying, we get
\begin{align*}
&D(x^T,p) - D(x,p)\\
&\leq \frac{4C_D}{T} + \frac{1}{(T-1)T} \sum_{t=1}^{T-3} t(t+1) \left(D(x^{t+1},p) - D(x^{t+1},p^t)\right)\\
&\quad + \frac{1}{(T-1)T} \left(\sum_{t=1}^{T-3} (t-1)t \left(D(x^t,p^t) - D(x^t,p)\right) + \sum_{t=1}^{T-3} \frac{4t}{t+1} \left(D(x,p^t)-D(x,p)\right) \right)\\
&\quad + \frac{T-2}{T} \left( \frac{T-3}{T-1} (D(x^{T-2},p^{T-2}) - D(x^{T-2},p)) + \frac{4}{(T-1)^2} (D(x,p^{T-2})-D(x,p)) \right)\\
&\quad + \frac{T-2}{T} (D(x^{T-1},p^{T-1}) - D(x^{T-1},p)) + \frac{4}{T^2} (D(x,p^{T-1})-D(x,p)).
\end{align*}
Applying Assumption \ref{ass:primal-distance-regular} completes the proof.
\Halmos
\endproof


\proof{Proof of Proposition~\ref{prop:infinite-curvature constant}.}

We will first show that the curvature constant $C_D$ defined in \eqref{eqn:curvature-constant} of $D(x,p) = \|x-p\|$ is infinite for any $p \in X$. Let us choose $x = p$, reserving the choice of $\alpha \in [0,1]$ and $s \in X$ for later. Then $D(x,p) = 0$, $D((1-\alpha)x + \alpha s,p) = \alpha \|s-p\|$, and the subgradients of $D(x,p)$ are $\left\{ y :~ \|y\|_* \leq 1 \right\}$. Thus, for any selection of subgradient mapping $y(\hat x) \in \grad_x D(\hat x,p)$ we have
\begin{align*}
\frac{1}{\alpha^2} \bigg[ D((1-\alpha)x + \alpha s,p) - D(x,p) - \alpha \la s-x, y(x) \ra \bigg] &= \frac{1}{\alpha^2} \bigg[ \alpha \|s-p\| - \alpha \la s-p, y(x) \ra \bigg]\\
&= \frac{1}{\alpha} \bigg[ \|s-p\| - \la s-p, y(x) \ra \bigg].
\end{align*}
Note that whenever there is a choice $s \in X$ with $\|s-p\| - \la s-p, y(p) \ra > 0$, we can send $\alpha \to 0$ and conclude that the curvature constant $C_D$ is infinite.

To choose the appropriate $s$, we denote the set of subgradients of $\|\cdot\|_q$ at $s-p$ as $G_{\|\cdot\|}(s-p)$. Observe that for a norm $\|\cdot\|$, if $y \in G_{\|\cdot\|}(s-p)$ then $\|y\|_* \leq 1$ and $\la s-p,y\ra = \|s-p\|$. Thus, we need to choose $s \in X$ such that $y(x) \not\in G_{\|\cdot\|}(s-p)$. To do this, we exploit the following property of $\ell_q$ norms. It is simple to check that for $q\in [1,\infty]$ and $y \in G_{\|\cdot\|_q}(s-p)$, we have the property that $y_{ij} > 0 \implies s_{ij} - p_{ij} > 0$. For our selection $y(x)$, first suppose that there exists $i \in \CA_j$ such that $y(x)_{ij} > 0$. Then a ranking $\sigma$ that ranks $i$ last will have $a(\sigma)_{ij} = 0$, so $a(\sigma)_{ij} - p_{ij} \leq 0$ because $p_{ij}\geq 0$. We cannot have $p = a(\sigma)$ for all $(n-1)!$ rankings $\sigma$ that ranks $i$ last (note that $n>2$); hence, there exists one $\sigma$ such that $a(\sigma) \neq p$, and we choose $s = a(\sigma)$. This implies that $y(x)_{ij} > 0$ while $s_{ij}-p_{ij} \leq 0$, hence $y(x) \not\in G_{\|\cdot\|_q}(s-p)$. Now suppose that $y(x)_{ij} \leq 0$ for all item-subset pairs $(i,j)$. If $y(x) = 0$, then the result follows trivially by choosing any $s \neq p$. Suppose now there exists some $y(x) < 0$. It is again simple to check that for $q\in [1,\infty]$ and $y \in G_{\|\cdot\|_q}(s-p)$, we have the property that $y_{ij} < 0 \implies s_{ij} - p_{ij} < 0$. Then a ranking $\sigma$ that ranks $i$ first will have $a(\sigma)_{ij} = 1$, so $a(\sigma)_{ij} - p_{ij} \geq 0$ because $p_{ij}\leq 1$. We cannot have $p = a(\sigma)$ for all $(n-1)!$ rankings $\sigma$ that ranks $i$ first, so there exists one such that $a(\sigma) \neq p$, and we choose $s = a(\sigma)$. This implies that $y(x)_{ij} < 0$ while $s_{ij}-p_{ij} \geq 0$, hence $y(x) \not\in G_{\|\cdot\|_q}(s-p)$. Thus, in all cases for $y(x)$, we can choose the appropriate $s \in X$.


In the case of MLE, we have $D(x,p) = -\sum_{j \in [m]} q_j \sum_{i \in \CA_j} p_{ij} \log(x_{ij}/p_{ij})$, which is the objective from \ref{sec:MLE-approach}. Here, $q_j$ is defined as in \eqref{eqn:pq_{ij}}. We can assume that $p_{ij} > 0$ by simply ignoring terms in the sum for which $p_{ij} = 0$. Choose $x=p$, which ensures that $D(\cdot,p)$ is differentiable at $x$ with $\grad_x D(x,p)_{ij} = -q_j /x_{ij}$. Then we have
\begin{align*}
&\frac{1}{\alpha^2}  \bigg[ D((1-\alpha)x + \alpha s,p) - D(x,p) - \alpha \la s-x, \grad_x D(x,p) \ra  \bigg]\\
&= -\frac{1}{\alpha^2} \sum_{j \in [m]} q_j \sum_{i \in \CA_j} p_{ij} \log\left(1-\alpha + \alpha \frac{s_{ij}}{p_{ij}}\right) + \frac{1}{\alpha} \sum_{j \in [m]} q_j \sum_{i \in \CA_j} \left(\frac{s_{ij}}{p_{ij}} - 1\right).
\end{align*}
Note that the second term is bounded by ${N\over \alpha} (\max_{i,j} 1/p_{ij} - 1)$. Choose $s_{ij} = a(\sigma)$ for any $\sigma \in S_n$. Then there exists some $i,j$ such that $s_{ij} = 0$. Sending $\alpha \to 1$ results in $\log\left(1-\alpha + \alpha \frac{s_{ij}}{p_{ij}}\right) \to \infty$, and the second term is bounded, so the curvature constant $C_D$ is infinite.
\Halmos
\endproof



\proof{Proof of Theorem~\ref{thm:online-choice-model-learning}.}
Note that we have
\begin{align*}
& D\left( x^T, p \right) = D\left( \frac{1}{T} \sum_{t=1}^{T} z^t,\, p \right) = \max_{y \in Y} \left\{ \left\la \frac{1}{T} \sum_{t=1}^{T} z^t, y \right\ra - D_*(y,p) \right\}\\
&= \max_{y \in Y} \frac{1}{T} \sum_{t=1}^T \left[ \la z^t, y \ra - D_*(y,p^t) + D_*(y,p^t) - D_*(y,p) \right]\\
&= \max_{y \in Y} \frac{1}{T} \sum_{t=1}^T \left[ f^t(y) + D_*(y,p^t) - D_*(y,p) \right]\\
&\leq \max_{y \in Y} \frac{1}{T} \sum_{t=1}^{T} f^t(y) - \frac{1}{T} \sum_{t=1}^T f^t(y^t) + \frac{1}{T} \sum_{t=1}^T f^t(y^t) + \max_{y \in Y} \frac{1}{T} \sum_{t=1}^T \left[ D_*(y,p^t) - D_*(y,p) \right],
\end{align*}
where the inequality follows from decomposing the terms in the maximum. Furthermore, we have $f^t(y^t) = \la z^t, y^t \ra - D_*(y^t;p^t) = f(y^t;p^t)$. Hence, from $y^t\in Y$, we deduce
\[ f^t(y^t) = f(y^t;p^t) \leq \max_{y \in Y} f(y,p^t) = \max_{y \in Y} \min_{x \in X} \left\{ \la x,y \ra - D_*(y,p^t) \right\} = \min_{x \in X} D(x,p^t). \]
Substituting this bound into the above expression and rearranging terms gives us
\begin{align*}
D\left( x^T, p \right) - \min_{x \in X} D(x,p) &\leq \max_{y \in Y} \frac{1}{T} \sum_{t=1}^{T} f^t(y) - \frac{1}{T} \sum_{t=1}^T f^t(y^t)\\
&\quad + \max_{y \in Y} \frac{1}{T} \sum_{t=1}^T \left[ D_*(y,p^t) - D_*(y,p) \right]\\
&\quad + \frac{1}{T} \sum_{t=1}^T \left[ \min_{x \in X} D(x,p^t) - \min_{x \in X} D(x,p) \right].
\end{align*}
Applying Assumption \ref{ass:dual-distance-regular} gives us the result.
\Halmos
\endproof

NonparChoice.bbl
\begin{thebibliography}{24}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\expandafter\ifx\csname urlstyle\endcsname\relax
  \providecommand{\doi}[1]{doi: #1}\else
  \providecommand{\doi}{doi: \begingroup \urlstyle{rm}\Url}\fi

\bibitem[Agrawal and Devanur(2015)]{AgrawalDevanur2015}
S.~Agrawal and N.~Devanur.
\newblock Fast algorithms for online stochastic convex programming.
\newblock In \emph{Proceedings of the Twenty-sixth Annual ACM-SIAM Symposium on
  Discrete Algorithms}, pages 1405--1424, Philadelphia, PA, USA, 2015. SIAM.

\bibitem[Ahmadi and Shanbhag(2014)]{AhmadiShanbhag2014}
H.~Ahmadi and U.~V. Shanbhag.
\newblock Data-driven first-order methods for misspecified convex optimization
  problems: Global convergence and rate estimates.
\newblock In \emph{53rd IEEE Conference on Decision and Control}, pages
  4228--4233, December 2014.

\bibitem[Bertsimas and Mi{\v{s}}ic(2015)]{BertsimasMisic2015}
D.~Bertsimas and V.~Mi{\v{s}}ic.
\newblock Data-driven assortment optimization.
\newblock Technical report, September 2015.
\newblock
  \url{https://www.algomus.com/wp-content/uploads/2016/11/Data-Driven-Assortment-Optimization.pdf}.

\bibitem[Block and Marschak(1960)]{BlockMarschak1960}
H.~Block and J.~Marschak.
\newblock Random orderings and stochastic theories of responses.
\newblock \emph{Contributions to Probability and Statistics}, 2:\penalty0
  97--132, 1960.

\bibitem[Calafiore and Campi(2005)]{CalafioreCampi2005}
G.~Calafiore and M.~Campi.
\newblock Uncertain convex programs: randomized solutions and confidence
  levels.
\newblock \emph{Mathematical Programming}, 102\penalty0 (1):\penalty0 25--46,
  2005.

\bibitem[Cesa-Bianchi and Lugosi(2006)]{CesaBianchiLugosi2006book}
N.~Cesa-Bianchi and G.~Lugosi.
\newblock \emph{Prediction, Learning, and Games}.
\newblock Cambridge University Press, New York, NY, USA, 2006.
\newblock ISBN 0521841089.

\bibitem[Desir et~al.(2016)Desir, Goyal, Jagabathula, and
  Segev]{DesirGoyalJagSegev2016}
A.~Desir, V.~Goyal, S.~Jagabathula, and D.~Segev.
\newblock Assortment optimization under the mallows model.
\newblock In D.~D. Lee, M.~Sugiyama, U.~V. Luxburg, I.~Guyon, and R.~Garnett,
  editors, \emph{Advances in Neural Information Processing Systems 29}, pages
  4700--4708. 2016.

\bibitem[Dwork et~al.(2001)Dwork, Kumar, Naor, and
  Sivakumar]{DworkKumarNaorSivakumar2001}
C.~Dwork, R.~Kumar, M.~Naor, and D.~Sivakumar.
\newblock Rank aggregation methods for the web.
\newblock In \emph{Proceedings of the 10th International Conference on World
  Wide Web}, WWW '01, pages 613--622, New York, NY, USA, 2001. ACM.

\bibitem[Farias et~al.(2009)Farias, Jagabathula, and Shah]{FariasJagShah2009}
V.~Farias, S.~Jagabathula, and D.~Shah.
\newblock A data-driven approach to modeling choice.
\newblock In Y.~Bengio, D.~Schuurmans, J.~D. Lafferty, C.~K.~I. Williams, and
  A.~Culotta, editors, \emph{Advances in Neural Information Processing Systems
  22}, pages 504--512. 2009.

\bibitem[Farias et~al.(2013)Farias, Jagabathula, and Shah]{FariasJS2013}
V.~Farias, S.~Jagabathula, and D.~Shah.
\newblock A nonparametric approach to modeling choice with limited data.
\newblock \emph{Management Science}, 59\penalty0 (2):\penalty0 305--322, 2013.

\bibitem[Hazan and Kale(2012)]{HazanKale2012}
E.~Hazan and S.~Kale.
\newblock Projection-free online learning.
\newblock In \emph{Proceedings of the 29th International Coference on
  International Conference on Machine Learning}, ICML'12, pages 1843--1850,
  USA, 2012. Omnipress.
\newblock ISBN 978-1-4503-1285-1.
\newblock URL \url{http://dl.acm.org/citation.cfm?id=3042573.3042808}.

\bibitem[Ho-Nguyen and K{\i}l{\i}n\c{c}-Karzan(2016)]{HoNguyenKK2016}
N.~Ho-Nguyen and F.~K{\i}l{\i}n\c{c}-Karzan.
\newblock Accelerating optimization under uncertainty via online convex
  optimization.
\newblock Technical report, August 2016.
\newblock \url{http://www.optimization-online.org/DB_HTML/2016/08/5571.html}.

\bibitem[Jagabathula and
  Rusmevichientong(2016)]{JagabathulaRusmevichientong2016}
S.~Jagabathula and P.~Rusmevichientong.
\newblock The limit of rationality in choice modeling: Formulation,
  computation, and implications.
\newblock Technical report, October 2016.
\newblock \url{https://ssrn.com/abstract=2859131}.

\bibitem[Jagabathula and Shah(2008)]{JagabathulaShah2008}
S.~Jagabathula and D.~Shah.
\newblock Inferring rankings under constrained sensing.
\newblock In D.~Koller, D.~Schuurmans, Y.~Bengio, and L.~Bottou, editors,
  \emph{Advances in Neural Information Processing Systems 21}, pages 753--760.
  2008.

\bibitem[Jaggi(2013)]{Jaggi2013}
M.~Jaggi.
\newblock Revisiting {Frank-Wolfe}: Projection-free sparse convex optimization.
\newblock In \emph{Proceedings of the 30th International Conference on Machine
  Learning}, volume~28 of \emph{Proceedings of Machine Learning Research},
  pages 427--435. PMLR, 2013.

\bibitem[Juditsky and Nemirovski(2012)]{JuditNem2012Pt1}
A.~Juditsky and A.~Nemirovski.
\newblock First-order methods for nonsmooth convex large-scale optimization,
  {I}: General purpose methods.
\newblock In S.~Sra, S.~Nowozin, and S.~Wright, editors, \emph{Optimization for
  Machine Learning}, Neural information processing series. MIT Press, 2012.

\bibitem[Mahajan and {van Ryzin}(2001)]{MahajanVanRyzin2001}
S.~Mahajan and G.~{van Ryzin}.
\newblock Stocking retail assortments under dynamic consumer substitution.
\newblock \emph{Operations Research}, 49\penalty0 (3):\penalty0 334--351, 2001.

\bibitem[Mirrokni et~al.(2017)Mirrokni, Leme, Vladu, and
  Wong]{MirrokniLemeVladuWong2017}
V.~Mirrokni, R.~P. Leme, A.~Vladu, and S.~C. Wong.
\newblock Tight bounds for approximate {C}arath{\'e}odory and beyond.
\newblock In \emph{Proceedings of the 34th International Conference on Machine
  Learning}, volume~70 of \emph{Proceedings of Machine Learning Research},
  pages 2440--2448, Aug 2017.

\bibitem[Nedi\'{c} and Ozdaglar(2009)]{NedicOzdaglar2009}
A.~Nedi\'{c} and A.~Ozdaglar.
\newblock Approximate primal solutions and rate analysis for dual subgradient
  methods.
\newblock \emph{SIAM Journal on Optimization}, 19\penalty0 (4):\penalty0
  1757--1780, 2009.

\bibitem[Nemirovski(2004)]{Nemirovski2004}
A.~Nemirovski.
\newblock Prox-method with rate of convergence {$O(1/t)$} for variational
  inequalities with {Lipschitz} continuous monotone operators and smooth
  convex-concave saddle point problems.
\newblock \emph{SIAM Journal on Optimization}, 15:\penalty0 229--251, 2004.

\bibitem[Rusmevichientong et~al.(2006)Rusmevichientong, Roy, and
  Glynn]{Rusmevichientong2006}
P.~Rusmevichientong, B.~V. Roy, and P.~Glynn.
\newblock A nonparametric approach to multiproduct pricing.
\newblock \emph{Operations Research}, 54\penalty0 (1):\penalty0 82--98, 2006.

\bibitem[Talluri and van Ryzin(2004)]{TallurivanRyzin2006book}
K.~Talluri and G.~van Ryzin.
\newblock \emph{The theory and practice of revenue management}, volume~68 of
  \emph{International Series in Operations Research \& Management Science}.
\newblock Springer US, 2004.

\bibitem[Train(2009)]{Train2009}
K.~E. Train.
\newblock \emph{Discrete Choice Methods with Simulation}.
\newblock Cambridge University Press, 2009.

\bibitem[van Ryzin and Vulcano(2015)]{vanRyzinVulcano2015}
G.~van Ryzin and G.~Vulcano.
\newblock A market discovery algorithm to estimate a general class of
  nonparametric choice models.
\newblock \emph{Management Science}, 61\penalty0 (2):\penalty0 281--300, 2015.

\end{thebibliography}
