\section{Introduction}

A multi-convex optimization problem is one in which the variables
can be partitioned into sets over which the problem is convex
when the other variables are fixed.
Multi-convex problems appear in domains such as
machine learning \cite{lee1999learning, udell2014generalized},
signal and information processing \cite{Lee00, Kim2008, Wen2012},
communication \cite{1254038},
and control \cite{751690, hassibi1999path, hours2014parametric, 7170815}.
Typical problems in these fields include nonnegative
matrix factorization (NMF) and bilinear matrix inequality (BMI) problems.

In general multi-convex problems are hard to solve globally,
but several algorithms have been proposed as heuristic or local methods,
and are widely used in applications.
Most of these methods are variations on
the block coordinate descent (BCD) method.
The idea of optimizing over a single block of variables while
holding the remaining
variables fixed in each iteration dates back to \cite{Warga,Powell1973}.
Convergence results were first discussed for strongly convex differentiable
objective function \cite{Warga}, and then under various assumptions
on the separability and regularity of the objective function
\cite{Tseng1993,Tseng2001}.
In \cite{attouch2010}, a two-block BCD method with proximal operator
is used to minimize a nonconvex objective function
which satisfies the Kurdyka-Lojasiewicz inequality.
In \cite{razaviyayn2012unified} the authors propose an inexact BCD approach which
updates variable blocks by minimizing a sequence of approximations
of the objective function, which can either be nondifferentiable or nonconvex.
A recent work \cite{xu2013block} uses BCD to solve multi-convex
problems, where the objective is a sum of a differentiable
multi-convex function and several extended-valued convex functions.
In each step, updates with and without proximal operator
and prox-linear operator are considered,
and convergence analysis is established under certain assumptions.
Gradient methods have also been proposed for multi-convex problems,
where the objective is differentiable in each block of variables,
and all variables are updated at once along their descent directions
and then projected into a convex feasible set in every iteration \cite{Chu04}.

The focus of this paper is not on solution methods,
but on a \emph{modeling framework} for expressing multi-convex
problems in a way that verifies the multi-convex structure, and
can expose the structure to whatever solution algorithm is then used.
Modeling frameworks have been developed for convex problems,
{\it e.g.},  \texttt{CVX} \cite{cvx}, \texttt{YALMIP} \cite{Lofberg:04},
\texttt{CVXPY} \cite{cvxpy_paper}, and \texttt{Convex.jl} \cite{cvxjl}.
These frameworks provide a uniform method for specifying convex
problems based on the idea of disciplined convex programming
(DCP) \cite{GBY:06}.
This gives a simple method for verifying that a problem is convex, and
for automatically canonicalizing to a standard
generic form such as a cone program.
The goal of DCP (and these software frameworks) is not to detect or
determine feasibility of an arbitrary problem, but rather to
give a very simple set of rules that can be used to construct
convex problems.

In this paper we extend the idea of DCP to multi-convex problems.
We propose a disciplined multi-convex programming (DMCP) rule set,
an extension of the DCP rule set.
Problem specifications that conform to the DMCP rule set can be verified as
convex in a group of variables, for any fixed values of the other variables,
using ideas that extend those in DCP.
We describe an efficient algorithm that can carry out the analysis of
convexity of problem when an arbitrary group of variables is fixed at any
value.
As with DCP, the goal of DMCP is not to analyze multi-convexity of
an arbitrary problem, but rather to give a simple set of rules which if followed
yields multi-convex problems.
In applications to date, such as NMF, verification of multi-convexity is simple, and
can be done by hand or just simple observation.
With DMCP a far larger class of multi-convex problems can be constructed
in an organized way.

We describe a software implementation of the ideas developed in this paper,
called \texttt{DMCP}, a Python package that extends \texttt{CVXPY}.
It implements the DMCP verification and analysis methods, and then
heuristically solves a conforming problem via BCD type algorithms,
which we extend for general use to include slack variables to handle
infeasibility.  A similar package, \texttt{MultiConvex}, has been
developed for the Julia package \texttt{Convex.jl}.
We illustrate the framework on a number of examples.

In \S\ref{s-multicvx} we carefully define multi-convexity of a function,
constraint, and problem.
We review block coordinate descent methods, introducing new variants
with slack variables and generalized inequalities, in \S\ref{s-bcd}.
In \S\ref{s-dmcp} we describe the main ideas of DMCP and an efficient
algorithm for verifying that a problem specification conforms to DMCP.
In \S\ref{s-impl} we describe our implementation of the package
\texttt{DMCP}.
Finally, in \S\ref{s-examples} we describe a number of numerical examples.
Our goal there is not to show competitive results, in terms of
solution quality or solve time, but rather to show the simplicity with which
the problem is specified, along with results that are at least
comparable to those obtained with custom solvers for the specific problem.

\section{Multi-convex programming}\label{s-multicvx}
\subsection{Multi-convex function}

\paragraph{Fixing variables in a function.}
Consider a function $f:{\mbox{\bf R}}^n\to{\mbox{\bf R}}\cup\{\infty\}$,
and a partition of the variable $x\in {\mbox{\bf R}}^n$ into blocks of variables
\[
x=(x_1,\ldots, x_N), \quad x_i\in{\mbox{\bf R}}^{n_i},
\]
so $\sum_{i=1}^N n_i = n$.
Throughout this paper we will use subsets of indices to refer to
sets of the variables.
Let $\mathcal{F}\subseteq\{1,\ldots,N\}$ denote an index set,
with complement $\mathcal{F}^c = \{1,\ldots,N\} \setminus \mathcal{F}$.
By fixing the variables with indices in $\mathcal{F}$
of the function $f$ at a given point $\hat{x}\in{\mbox{\bf R}}^n$,
we obtain a function over the remaining variables, with indices in $\mathcal{F}^c$,
which we denote as $\tilde f = \mathrm{fix}(f, \hat{x}, \mathcal{F}) $.
For $i \in \mathcal{F}^c$, $x_i$ is a variable of the function $\tilde f$;
for $i \in \mathcal{F}$, $x_i=\hat x_i$.
Informally we refer to $\tilde f$ as
`$f$, with the variables $x_i$ for $i\in \mathcal{F}$ fixed'.

As an example consider $f:{\mbox{\bf R}}^4 \to {\mbox{\bf R}} \cup \{\infty\}$ defined as
\begin{equation}\label{e-example}
f(x_1, x_2, x_3,x_4) = |x_1x_2 + x_3x_4|.
\end{equation}
With $\hat{x} = (1,2,1,3)$ and $\mathcal{F}=\{1,3\}$,
the fixed function $\tilde f =
\mathrm{fix}(f, \hat{x}, \{1,3\})$ is given by
$\tilde f(x_2, x_4)  = |x_2+x_4|$.

\paragraph{Multi-convex and multi-affine functions.}
Given an index set $\mathcal{F}\subseteq\{1,\ldots,N\}$,
we say that a function $f:{\mbox{\bf R}}^n\to{\mbox{\bf R}}\cup\{\infty\}$
is \textit{convex} (or \textit{affine}) \textit{with set $\mathcal{F}$ fixed},
if for any $\hat{x}\in{\mbox{\bf R}}^n$
the function $\mathrm{fix}(f,\hat{x},\mathcal{F})$ is a convex (or affine) function.
(In this definition, we consider a so-called improper function
\cite{rockafellar}, which has the value $\infty$ everywhere,
as convex or affine.)
For example, the function $f$ defined in \eqref{e-example}
is convex with the variables $x_1$ and $x_3$ fixed ({\it i.e.},
with index set $\mathcal{F}=\{1,3\}$).
A function is convex if and only if it is convex
with $\mathcal{F} = \emptyset$ fixed,
{\it i.e.}, with none of its variables fixed.

We say the function $f$ is \textit{multi-convex}
(or \textit{multi-affine)}, if there are index sets $\mathcal{F}_1,\ldots,\mathcal{F}_K$,
such that for every $k$ the function $f$ is convex (or affine) with $\mathcal{F}_k$ fixed,
and $\cap_{k=1}^K \mathcal{F}_k = \emptyset$.
The requirement that $\cap_{k=1}^K \mathcal{F}_k = \emptyset$ means that for
every variable $x_i$ there is some $\mathcal{F}_k$ with $i \not\in \mathcal{F}_k$.
In particular, $\tilde{f}= \mathrm{fix}(f,\hat{x},\{1,\ldots,N\}\setminus \{i\})$
is convex in $x_i$.
For $K=2$, we say that the function is \emph{bi-convex}
(or \emph{bi-affine}).

As an example, the function $f$ in \eqref{e-example} is
convex with $\{1,3\}$ fixed and $\{2,4\}$ fixed, so it is multi-convex.
The choice of $K$, and the index sets, is not unique.
The function $f$ is also convex
with $\{2,3,4\}$ fixed, $\{1,3,4\}$ fixed, $\{1,2,4\}$ fixed, and $\{1,2,3\}$ fixed.

\paragraph{Minimal fixed sets.}
For a function $f$, we can consider the set of all index sets $\mathcal{F}$ for which
$\mathrm{fix}(f, \hat x, \mathcal{F})$ is convex for all $\hat x$;
among these we are interested in
the \emph{minimal fixed sets} that render a function convex.
A minimal fixed set is a set of variables that when fixed make the
function convex; but if any variable is removed from the set, the function
is not convex.
A function is multi-convex if and only if the intersection of these
minimal fixed index sets is empty.

\subsection{Multi-convex problem}
We now extend the idea of multi-convexity to the optimization problem
\begin{equation}\label{main}
\begin{array}{ll}
\mbox{minimize} &  f_0(x) \\
\mbox{subject to} & f_i (x) \leq 0, \quad i = 1,\ldots, m \\
& g_i (x)  = 0, \quad i = 1,\ldots, p,
\end{array}
\end{equation}
with variable $x\in{\mbox{\bf R}}^n$ partitioned into blocks as
$x = (x_1,\ldots,x_N)$,
and functions $f_i:{\mbox{\bf R}}^{n}\to{\mbox{\bf R}}\cup\{\infty\}$ for $i=0,\ldots,m$
and $g_i:{\mbox{\bf R}}^n\to{\mbox{\bf R}}$ for $i=1,\ldots, p$ are proper.

Given an index set $\mathcal{F}\subseteq\{1,\ldots,N\}$,
problem~\eqref{main} is \textit{convex with set $\mathcal{F}$ fixed},
if for any $\hat{x}\in{\mbox{\bf R}}^n$ the problem
\begin{equation}\label{fix_p}
\begin{array}{ll}
\mbox{minimize} &  \mathrm{fix}(f_0, \hat{x}, \mathcal{F}) \\
\mbox{subject to} &\mathrm{fix} (f_i, \hat{x}, \mathcal{F}) \leq 0, \quad i = 1,\ldots, m \\
& \mathrm{fix} (g_i, \hat{x}, \mathcal{F})   = 0, \quad i = 1,\ldots, p,
\end{array}
\end{equation}
is convex.
In other words, problem~\eqref{main} is convex with $\mathcal{F}$ fixed,
if and only if functions $f_i$ for $i=0,\ldots, m$ are convex with $\mathcal{F}$ fixed,
and functions $g_i$ for $i=1,\ldots,p$ are affine with $\mathcal{F}$ fixed.

We say the problem~\eqref{main} is \textit{multi-convex}, if there are sets
$\mathcal{F}_1,\ldots, \mathcal{F}_K$,
such that for every $k$ problem~\eqref{main} is convex with set $\mathcal{F}_k$ fixed,
and $\cap_{k=1}^K \mathcal{F}_k = \emptyset$.
A convex problem is multi-convex with $K=0$ ({\it i.e.}, $\mathcal{F} = \emptyset$).
A \emph{bi-convex problem} is multi-convex with $K=2$.

As an example the following problem is multi-convex:
\begin{equation}\label{example}
\begin{array}{ll}
\mbox{minimize} & |x_1x_2 + x_3x_4| \\
\mbox{subject to} & x_1+x_2+x_3+x_4 = 1,
\end{array}
\end{equation}
with variable $x\in {\mbox{\bf R}}^4$.
This is readily verfied with
$\mathcal F_1 = \{1,3\}$ and $\mathcal F_2 = \{2,4\}$.
For a given problem we can consider the minimal variable index
sets which make the problem convex.  If the problem is convex,
$\mathcal{F} = \emptyset$ is the unique minimal set.

\section{Block coordinate descent and variations}\label{s-bcd}

In this section we review, and extend, some generic methods for approximately
solving the multi-convex problem \eqref{main}, using BCD-type methods.

\subsection{Block coordinate minimization with slack variables}
Assume that sets $\mathcal{F}_k$, $k=1,\ldots,K$ are index sets
for which the problem~\eqref{main} with $\mathcal{F}_k$ fixed is convex,
with $\cap_{k=1}^K \mathcal F_k =\emptyset$.
These could be the set of all minimal index sets, but any other set of index sets
that verify multi-convexity could be used.

The basic form of the proposed method is iterative.  In each iteration,
we fix the variables in one set $\mathcal{F}_k$ and solve the following subproblem,
\begin{equation}\label{slack}
\begin{array}{ll}
\mbox{minimize} &  \mathrm{fix}( f_0, \hat{x}, \mathcal{F}_k)
+\mu \sum_{i=1}^{m} s_i + \mu\sum_{i=1}^p |s_{i+m}|\\
\mbox{subject to} &  \mathrm{fix}( f_i, \hat{x}, \mathcal{F}_k) \leq s_i, \quad s_i\geq0,\quad i = 1,\ldots, m \\
&  \mathrm{fix}( g_i, \hat{x}, \mathcal{F}_k)  = s_{i+m}, \quad i = 1,\ldots, p, \\
\end{array}
\end{equation}
where $s_i$ for $i=1,\ldots,m+p$ and $x_i$ for $i\in\mathcal{F}_k^c$ are the variables,
and $\mu>0$ is a parameter.
Here the constant $\hat{x}$ inherits the value of $x$ from the least iteration.
This subproblem solved in each iteration is convex.
The slack variables $s_i$ for $i=1,\ldots,m+p$ ensure that the subproblem
\eqref{slack} cannot be infeasible.
The added terms in the objective are a so-called exact penalty
\cite{nocedal2006numerical},
meaning that when some technical conditions hold, and $\mu$ is large enough,
the solution satisfies $s_i=0$, when the subproblem without the slack
variables is feasible.

Many schemes can be used to choose $k$ in each iteration, and to update the
slack parameter $\mu$.
For example, we can cyclically choose $k$, or randomly choose $k$,
or optimize over $k=1, \ldots, K$ in rounds of $K$ steps,
in an order chosen by a random permutation in each round.
Updating $\mu$ is typically done by increasing it by a factor $\rho>1$ after each
iteration, or after each round of $K$ iterations.
One variation on the algorithm sets the slack variables to zero ({\it i.e.}, removes
them) once a feasible point is obtained ({\it i.e.}, a point is obtained with
$s_i=0$).
The algorithm is typically initialized with values specific to the particular
application, or generic values.

This algorithm differs from the basic BCD algorithm in the addition
of the slack variables, and in the feature that a variable can appear in more
than one set $\mathcal{F}_k^c$, meaning that a variable can be updated
in multiple iterations per round of $K$ iterations.
For example, if a problem has variables $x_1,x_2,x_3$ and is convex in
$(x_1,x_2)$ and $(x_2,x_3)$, our method will update $x_2$ in each step.

In the general case, very little can be said about the convergence of this method.
One obvious observation is that, if $\mu$ is held fixed, the objective
is nonincreasing and so convergences.
See the references cited above for some convergence results for related algorithms,
for special cases with strong assumptions such as strict convexity (when the variables
are fixed) or differentiability.
As a practical matter, similar algorithms have been found to be robust, and
very useful in practice, despite a lack of strong theory
establishing convergence in the general case.

\subsection{Block coordinate proximal iteration}
A variation of subproblem \eqref{slack} is adds a proximal term
\cite{parikh2014proximal},
which renders the subproblems strongly convex:
\begin{equation}\label{proximal}
\begin{array}{ll}
\mbox{minimize} &   \mathrm{fix}( f_0, \hat{x}, \mathcal{F}_k)
+\mu \sum_{i=1}^{m} s_i + \mu \sum_{i=1}^p |s_{i+m}|
+\frac{1}{2\lambda} \sum_{i\in\mathcal{F}^c_k} \|x_i-\hat{x}_i\|_2^2\\
\mbox{subject to}
& \mathrm{fix}( f_i, \hat{x}, \mathcal{F}_k) \leq s_i, \quad s_i\geq0, \quad i = 1,\ldots, m \\
&  \mathrm{fix}( g_i, \hat{x}, \mathcal{F}_k) = s_{i+m}, \quad i = 1,\ldots, p,
\end{array}
\end{equation}
where $x_i$ for $i\in\mathcal{F}^c_k$ and $s_i$ for $i=1,\ldots, m+p$ are variables,
$\lambda >0$ is the proximal parameter.
The proximal term penalizes large changes in the variables being
optimized, {\it i.e.}, it introduces damping into the algorithm.  In some cases
it has been observed to yield better final points, {\it i.e.}, points
with smaller objective value, than those obtained without proximal
regularization.

Yet another variation uses linearized proximal steps, when $f$ is differentiable
in the variables $x_i$ for $i \in \mathcal F_k^c$.  The subproblem
solved in this case is
\begin{equation}\label{linear}
\begin{array}{ll}
\mbox{minimize} &
\mu \sum_{i=1}^{m} s_i + \mu\sum_{i=1}^p |s_{i+m}|
+ \sum_{i\in\mathcal{F}^c_k} \left(  \frac{1}{2\lambda} \|x_i-\hat{x}_i\|_2^2
+(x_i-\hat{x}_i)^T \nabla f(\hat{x}_i) \right) \\
\mbox{subject to}
&  \mathrm{fix}( f_i, \hat{x}, \mathcal{F}_k)\leq s_i, \quad s_i\geq0, \quad i = 1,\ldots, m \\
& \mathrm{fix}( g_i, \hat{x}, \mathcal{F}_k)  = s_{i+m}, \quad i = 1,\ldots, p,
\end{array}
\end{equation}
where $x_i$ for $i\in\mathcal{F}^c_k$ and $s_i$ for $i=1,\ldots, m+p$ are variables,
and $\nabla f(\hat{x}_i)$ is the partial gradient of $f$ with respect to $x_i$ at
the point $\hat{x}$.
The objective is equivalent to the minimization of
\[
\mu \sum_{i=1}^{m} s_i + \mu\sum_{i=1}^p |s_{i+m}|
+\sum_{i\in\mathcal{F}^c_k} \frac{1}{2\lambda} \|x_i-\hat{x}_i +
\lambda \nabla f(\hat{x}_i)\|_2^2,
\]
which is the objective of a proximal gradient method.

\subsection{Generalized inequality constraints}
One useful extension is to generalize problem~\eqref{main} by allowing
generalized inequality constraints.
Suppose the functions $f_0$ and $g_i$ for $i=1,\ldots, p$ are the same as
in problem~\eqref{main},
but $f_i : {\mbox{\bf R}}^n\to{\mbox{\bf R}}^{d_i}\cup\{\infty\}$.
Consider the following program with generalized inequalities,
\begin{equation}
\begin{array}{ll}
\mbox{minimize} & f_0(x) \\
\mbox{subject to} & f_i (x) \preceq_{\mathcal{K}_i} 0,\quad i=1,\ldots, m \\
& g_i(x) = 0,\quad i = 1,\ldots,p,
\end{array}
\end{equation}
where $x = (x_1,\ldots, x_N)\in{\mbox{\bf R}}^n$ is the variable,
and the generalized inequality constraints
are with respect to proper cones $\mathcal{K}_i\subseteq{\mbox{\bf R}}^{d_i}$,
$i=1, \ldots, m$.
The definitions of multi-convex program and minimal index set
can be directly extended.
Slack variables are added in the following way:
\begin{equation}\label{ext}
\begin{array}{ll}
\mbox{minimize} &  \mathrm{fix}( f_0, \hat{x}, \mathcal{F}_k)
+\mu \sum_{i=1}^{m} s_i + \mu\sum_{i=1}^p |s_{i+m}|\\
\mbox{subject to}
&  \mathrm{fix}( f_i, \hat{x}, \mathcal{F}_k) \preceq_{\mathcal{K}_i} s_i e_i,
 \quad s_i\geq0,\quad i = 1,\ldots, m \\
&  \mathrm{fix}( g_i, \hat{x}, \mathcal{F}_k)  = s_{i+m}, \quad  i = 1,\ldots, p,
\end{array}
\end{equation}
where $e_i$ is a given positive element in cone $\mathcal{K}_i$
for $i=1,\ldots,m$.

\section{Disciplined multi-convex programming}\label{s-dmcp}
\subsection{Disciplined convex programming}
Disciplined convex programming (DCP) is a methodology introduced by Grant et al.
\cite{GBY:06}
that imposes a set of conventions that must be followed when
constructing (or specifying or defining) convex programs.
Conforming problems are called \emph{disciplined convex programs}.
A disciplined convex program can be transformed into an
equivalent cone program by replacing each function
with its graph implementation \cite{GB:08}.
The convex optimization modeling systems
\texttt{YALMIP} \cite{Lofberg:04}, \texttt{CVX} \cite{cvx},
\texttt{CVXPY} \cite{cvxpy_paper},
and \texttt{Convex.jl} \cite{cvxjl} use DCP
to verify the convexity of a problem
and automatically convert convex programs into cone programs,
which can then be solved using generic solvers.

The conventions of DCP restrict the set of functions that can appear in
a problem and the way functions can be composed.
Every function in a disciplined convex program must be formed
as an expression involving constants or parameters, variables,
and a dictionary of atomic functions.
The dictionary consists of functions with known curvature and monotonicity,
and a graph implementation, or representation as partial
optimization over a cone program \cite{BoV:04,NesNem:92}.
Every composition of functions $f(g_1(x),\ldots,g_p(x))$,
where $f : {\mbox{\bf R}}^p \to {\mbox{\bf R}}$ is convex and
$g_1,\ldots,g_p : {\mbox{\bf R}}^n \to {\mbox{\bf R}}$,
must satisfy the following composition rule,
which ensures the composition is convex.
Let $\tilde{f} : {\mbox{\bf R}}^p \to {\mbox{\bf R}} \cup \{\infty\}$
be the extended-value extension of $f$ \cite[Chap.~3]{BoV:04}.
One of the following conditions must hold for each $i=1,\ldots,p$:
\begin{itemize}
	\item $g_i$ is convex and $\tilde{f}$ is nondecreasing in argument $i$.
	\item $g_i$ is concave and $\tilde{f}$ is nonincreasing in argument $i$.
	\item $g_i$ is affine.
\end{itemize}
The composition rule for concave functions is analogous.

Signed DCP is an extension of DCP that keeps track of the
signs of functions and expressions, using simple sign arithmetic.
The monotonicity of functions in the atom library can then depend
on the sign of their arguments.
As a simple example,
consider the expression $y=(\exp x)^2$, where $x$ is a variable.
The subexpression $\exp x$ is convex and (in signed DCP analysis)
nonnegative.  With DCP analysis, $y$ cannot be verified as
convex, since the square function is not nondecreasing. With signed
DCP analysis, the square function is known to be nondecreasing
for nonnegative arguments, which matches this case, so $y$ is
verified as convex using signed DCP analysis.

Convexity verification of an expression formed from variables
and constants (or parameters) in (signed) DCP first analyzes the signs
of all subexpressions.  Then it checks that the composition rule
above holds for every subexpression, possibly relying on the known
signs of subexpressions.
If everything checks out the expression is verified to be constant,
affine, convex, concave, or unknown (when the DCP rules do not hold).
We make an observation that is critical for our work here:
The DCP analysis does not use the values of any constants or parameters
in the expression.  The number 4.57 is simply treated as positive;
if a parameter has been declared as positive, then it is treated
as positive.
It follows immediately that DCP analysis has verified not just that
the specific expression is convex, but that it is convex for
\emph{any} other values of the constants and parameters, with the
same signs as the given ones, if the sign matters.

The following code snippet gives an example in \texttt{CVXPY}:
\begin{quote}
	\begin{verbatim}
	x = Variable(n)
	mu = Parameter(sign = 'positive')
	expr = sum_squares(x) + mu*norm(x,1)
	\end{verbatim}
\end{quote}
In the first line we declare (or construct) a variable,
and in the second line construct a paramater, {\it i.e.}, a constant
that is unknown, but positive.
The curvature of the expression \texttt{expr} is verified to be convex,
even though the value of parameter \texttt{mu} is unknown;
DCP analysis uses only the fact that whatever the value of $\texttt{mu}$
is, it must be positive.

\subsection{Disciplined multi-convex programming}

To determine that problem~\eqref{main}
is multi-convex requires us to verify that functions
$\mathrm{fix}(f_i,\hat{x},\mathcal{F})$ are convex,
and if $\mathrm{fix}(g_i,\hat{x},\mathcal{F})$ are affine, for all $\hat{x}\in{\mbox{\bf R}}^n$.
We can use the idea of DCP, specifically with signed parameters,
to carry this out,
which gives us a practical method for multi-convexity verification.

\paragraph{Multi-convex atoms.}
We start by generalizing the library of DCP atom functions to include
multi-convex atomic functions.
A function is a multi-convex atom if it has $N$ arguments $N>1$,
and it reduces to a DCP atomic function, when and only when
all but the $i$th arguments are constant for each $i=1,\ldots,N$.
For example, the product of $N$ variables is a multi-convex atom
that extends the DCP atom of multiplication between one variable and constants.

Given a description of problem~\eqref{main} under a library of DCP and
multi-convex atomic functions,
we say that it is \textit{disciplined convex programming
	with set $\mathcal{F}\subseteq\{1,\ldots,N\}$ fixed},
if the corresponding problem~\eqref{fix_p} for any $\hat{x}\in{\mbox{\bf R}}^n$
conforms to the DCP rules with respect to the DCP atomic function set.
When there is no confusion,
we simply say that problem~\eqref{main} is DCP with $\mathcal{F}$ fixed.

To verify if a problem is DCP with $\mathcal{F}$ fixed,
a method first fixes the problem by replacing variables in $\mathcal{F}$
with parameters of the same signs and dimensions.
Then it verifies DCP of the fixed problem with parameters
according to the DCP ruleset.
The parameter is the correct model for fixed variables,
in that the DCP rules ensure that the verified curvature holds for
any value of the parameter.

\paragraph{Disciplined multi-convex program.}
Given a description of problem~\eqref{main} under a library of DCP and
multi-convex atomic functions,
it is \textit{disciplined multi-convex programming} (DMCP),
if there are sets $\mathcal{F}_1,\ldots, \mathcal{F}_K$ such that
problem~\eqref{main} with every $\mathcal{F}_k$ fixed is DCP,
$\cap_{k=1}^K \mathcal{F}_k = \emptyset$.
We simply say that problem~\eqref{main} is DMCP if there is no confusion.
A problem that is DMCP is guaranteed to be multi-convex, just as
a problem that is DCP is guaranteed to be convex.
Morever, when a BCD method is applied to a DMCP problem, each iteration
involves the solution of a DCP problem.

\paragraph{DMCP verification.}
A direct way of DMCP verification is to check if the problem is DCP with
$\{i\}^c$ fixed for every $i=1,\ldots,N$.
Expressions in DMCP inherit the tree structure from DCP,
so such verification can be done in $O(MN)$ time,
where $M$ is number of nodes in problem expression trees,
and $N$ is the number of distinct variables.
To see why DMCP can be verified in this simple way,
we have the following claim.
The claim implies that every DMCP problem is DCP when all but one variable
is fixed.

\begin{claim}\label{claim:DCP_exp}
	For a problem consisting only of DCP atoms (or multi-convex atoms
	with all but one arguments constant),
	it is DCP with $\mathcal{F}$ fixed,
	if and only if it is DCP with $\{i\}^c$ fixed for all $i\in\mathcal{F}^c$.
\end{claim}

To see why Claim~\ref{claim:DCP_exp} is correct, we first prove the direction
that if a problem consisting only of DCP atoms is DCP with $\{i\}^c$ fixed
for all $i\in\mathcal{F}^c$, then it is DCP with $\mathcal{F}$ fixed.
The proof begins with two observations for functions consisting only of DCP
atoms.
\begin{itemize}
	\item The DCP curvature types have a hierarchy:
	\[
	\mathrm{unknown}\to\left\{
	\begin{array}{ll}
	&\mathrm{convex} \\
	&\mathrm{concave}
	\end{array}\;\;
	\right\} \to \mathrm{affine} \to \mathrm{constant}.
	\]
	Unknown is the base type, then it splits into convex and concave.
	Affine is a subtype of both convex and concave.
	Then constant is a subtype of everything.
	There is a similar hierarchy for sign information.
	The DCP type system is monotone in the curvature and sign hierarchies,
	meaning if the curvature or sign of an argument of a function is changed
	to be more specific, the type of the function will become more specific
	or stay the same.
	Fixing variables of a function
	makes the curvatures of some arguments more specific,
	while keeps all signs the same,
	so the curvature of the function can only get more specific or stay the same.

	\item If a function is affine in $x$ and $y$ separately,
	then it is affine in $(x,y)$.
	This is true because no multiplication of variables is allowed in DCP,
	and the function can only be in the form of $Ax+By+c$.
\end{itemize}
Now suppose that a problem consisting of DCP atoms
with only two variables $x$ and $y$ is not DCP,
but it is DCP with $x$ fixed and with $y$ fixed.
Then there must be some function that has a wrong curvature type in $(x,y)$
but whose arguments all have known curvatures.
Since the function has the right curvature type with $x$ fixed and $y$ fixed,
there must be an argument that is convex or concave (not affine) in $(x,y)$,
but has a different curvature in $x$ and $y$.
According to the first observation, the argument must be affine in $x$ and $y$.
By the second observation, the argument is affine in $(x,y)$,
which is a contradiction.
For the same reason, cases with more than two variables have the same conclusion.

For the other direction of Claim~\ref{claim:DCP_exp},
we again use the observation that if a problem is DCP with $\mathcal{F}$ fixed,
then fixing additional variables only makes function curvatures
more specific.
The problem must then be DCP with $\{i\}^c$ fixed
for all $i \in \mathcal{F}^c$, since $\mathcal{F} \subseteq \{i\}^c$.

\subsection{Efficient search for minimal sets}\label{s-min-sets}
A problem may have multiple collections of index sets for which it is DMCP.
We propose several generic and efficient ways of choosing which collection to use
when applying a BCD method to the problem.
The simplest option is to always choose the collection $\{1\}^c,\ldots,\{N\}^c$,
in which case BCD optimizes over one variable at a time.

A more sophisticated approach is to reduce the collection $\{1\}^c,\ldots,\{N\}^c$ to
minimal sets,
which allows BCD to optimize over multiple variables each iteration.
We find minimal sets by first determining which variables can be optimized
together.
Concretely, we construct a conflict graph $(\mathcal{V}, \mathcal{E})$,
where $\mathcal{V}$ is the set of all variables,
and $i\sim j\in\mathcal{E}$ if and only if variables $i$ and $j$
appear in two different child trees of a multi-convex atom in the problem
expression tree,
which means the variables cannot be optimized together.

Constructing the conflict graph takes $O(N^2M)$ time.
We simply do a depth-first traversal of the problem expression tree.
At each leaf node, we initialize a linked list with the leaf variable.
At each parent node, we join the linked lists of its children.
At each multi-convex atom node, we also remove duplicates from each child's
linked list and then iterate over the lists,
adding an edge for every two variables appearing in different lists.
The edges added at a given multi-convex node are all unique
because a duplicate edge would mean the same variable appeared in two different
child trees, which is not possible in a DMCP problem.
Hence, iterating over the lists of variables takes at most $N^2$ operations.

Given the conflict graph, for $i=1,\ldots,N$ we find a maximal independent
set $\mathcal{F}_i$ containing $i$ using a standard fast algorithm and replace $\{i\}^c$
with $\mathcal{F}_i^c$.
The final collection is all index sets $\mathcal{F}_i^c$ that are not supersets
of another index set $\mathcal{F}_j^c$.
More generally, we can choose any collection of index sets $\mathcal{F}_1,\ldots,\mathcal{F}_K$ such
that $\mathcal{F}_1^c,\ldots,\mathcal{F}_K^c$ are independent sets in the conflict graph.

\section{Implementation}\label{s-impl}
The methods of DMCP verification, searching for minimal sets to fix,
and cyclic optimization with
minimal sets fixed are implemented as an extension of \texttt{CVXPY}
in a package \texttt{DMCP} that can be accessed at
\verb|https://github.com/cvxgrp/dmcp|.
A Julia package with similar functionality, \texttt{MultiConvex.jl},
can be found at
\verb|https://github.com/madeleineudell/MultiConvex.jl|,
but we focus here on the Python package \texttt{DMCP}.

\subsection{Some useful functions}
\paragraph{Multi-convex atomic functions.}
In order to allow multi-convex functions,
we extend the atomic function set of \texttt{CVXPY}.
The following atoms are allowed to have non-constant expressions
in both arguments,
while in base \texttt{CVXPY} one of the arguments must be constant.
\begin{itemize}
\item
multiplication: \verb|expression1 * expression2|
\item
elementwise multiplication: \verb|mul_elemwise(expression1, expression2)|
\item
convolution: \verb|conv(expression1, expression2)|
\end{itemize}

\paragraph{Find minimal sets.}
Given a problem, the function
\begin{quote}
	\begin{verbatim}
	find_minimal_sets(problem)
	\end{verbatim}
\end{quote}
runs the algorithm discussed in \S~\ref{s-min-sets}
and returns a list of minimal sets of indices of variables.
The indices are with respect to the list \verb|problem.variables()|,
namely, the variable corresponding to index $0$ is
\verb|problem.variables()[0]|.

\paragraph{DMCP verification.}
Given a problem, the function
\begin{quote}
	\begin{verbatim}
	is_dmcp(problem)
	\end{verbatim}
\end{quote}
returns a boolean indicating if it is a DMCP problem.

\paragraph{Fix variables.}
The function
\begin{quote}
	\begin{verbatim}
	fix(expression, fix_vars)
	\end{verbatim}
\end{quote}
returns a new expression with the variables in the list \verb|fix_vars|
replaced with parameters of the same signs and values.
If \verb|expression| is replaced with a \verb|CVXPY| problem,
 then a fixed problem is returned by fixing every expression in its
 cost function and both sides of inequalities or equalities in constraints.

\paragraph{Random initialization.}
It is suggested that users provide an initial point $x^0$ for the method such that
functions $\mathrm{fix}(f_i,x^0,\mathcal{F}_1)$ are proper for $i=0,\ldots,m$,
where $\mathcal{F}_1$ is the first minimal set given by \verb|find_minimal_sets|.
If not, the function \verb|rand_initial(problem)|
will be called to generate random values from the uniform distribution
over the interval $[0,1)$ ($(-1,0]$) for variables with non-negative (non-positive) sign,
and from the standard normal distribution for variables with no sign.
There is no guarantee
that such a simple random initialization can always work for any problem.

\subsection{Options on update and algorithm parameters}
The solving method is to cyclically fix every minimal set
found by \verb|find_minimal_sets| and update the variables.
Three ways of updating variables are implemented.
The default one can be called by
\verb|problem.solve(method = 'bcd', update = 'proximal')|,
which is to solve the subproblem with proximal operators,
{\it i.e.}, problem~\eqref{proximal}.
To update by minimizing the subproblem without proximal operators,
{\it i.e.}, problem~\eqref{slack}, the solve method is called with
\verb|update = 'minimize'|.
To use the prox-linear operator in updates, {\it i.e.}, problem~\eqref{linear},
the solve method should be called with
\verb|update = 'prox_linear'|.

The parameter $\mu$ is update in every cycle by
$\mu_{t+1} = \mathrm{min} (\rho\mu_t, \mu_\mathrm{max})$.
The algorithm parameters are $\rho,\mu_0, \mu_\mathrm{max}, \lambda$,
and the maximum number of iterations.
They can be set by passing values of the parameters
 \verb|rho|, \verb|mu_0|, \verb|mu_max|,
\verb|lambd|, and \verb|max_iter|, respectively, to the \verb|solve| method.

\section{Numerical examples}\label{s-examples}
\subsection{One basic example}
\paragraph{Problem description.}
The first example is problem~\eqref{example},
which has appeared throughout this paper to explain definitions.
\paragraph{DMCP specification.}
The code written in \verb|DMCP| for this example is as follows.
\begin{quote}
	\begin{verbatim}
	prob = Problem(Minimize(abs(x_1*x_2+x_3*x_4)), [x_1+x_2+x_3+x_4 == 1])
	\end{verbatim}
\end{quote}
To find all minimal sets, the following line is typed in
\begin{quote}
	\begin{verbatim}
	find_minimal_sets(prob)
	\end{verbatim}
\end{quote}
and the output is \texttt{[[2, 1], [3, 1], [2, 0], [3, 0]]}.
Note that index $i$ corresponds to variable \verb|prob.variables()[i]|
for $i=0,\ldots, 3$.
To verify if it is DMCP, the function
\begin{quote}
	\begin{verbatim}
	is_dmcp(prob)
	\end{verbatim}
\end{quote}
returns \verb|True|.

\paragraph{Numerical result.}

Random initial values are set for all variables. The solve method with
default setting  finds a feasible point
with objective value $0$, which solves the problem globally.

\subsection{Fractional optimization}

\paragraph{Problem description.}
In this example we evaluate DMCP on some fractional optimization problems.
Consider the following problem
\begin{equation}\label{frac}
\begin{array}{ll}
\mbox{minimize} & p(x)/q(x) \\
\mbox{subject to} & x\in\mathcal{X},
\end{array}
\end{equation}
where $x\in{\mbox{\bf R}}^n$ is the variable, $\mathcal{X}$ is a convex set,
$p$ is a convex function, and $q$ is concave. The objective function
is set to $+\infty$ unless $p(x) \geq 0, q(x)>0$.
Such a problem is quasi-convex, and can be globally solved \cite[\S4.2.5]{BoV:04},
or even have analytical solutions,
so the aim here is just to evaluate the effectiveness of the method.

There are several ways of specifying problem~\eqref{frac} as DMCP.
One way is via the following problem.
\begin{equation}\label{frac_t1}
\begin{array}{ll}
\mbox{minimize} & p(x)/q(y) \\
\mbox{subject to} & x\in\mathcal{X}, \quad x = y,
\end{array}
\end{equation}
where $x$ and $y$ are variables.
Another way is via the following.
\begin{equation}\label{frac_t2}
\begin{array}{ll}
\mbox{minimize} & \alpha \\
\mbox{subject to} & x\in\mathcal{X}, \quad p(x) \leq \alpha q(x),
\end{array}
\end{equation}
where $\alpha\in{\mbox{\bf R}}_+$ and $x$ are variables.
Both of them are biconvex.

\paragraph{DMCP specification.}
Suppose that $\mathcal{X} = {\mbox{\bf R}}^n$.
The code for the formulation in \eqref{frac_t1} is as the following.
\begin{quote}
	\begin{verbatim}
	x = Variable(n)
	y = Variable(n)
	# specify p and q here
	prob = Problem(Minimize(inv_pos(q)*p), [x == y])
	\end{verbatim}
\end{quote}
Expressions \texttt{p} and \texttt{q} are to be specified.
The code for problem~\eqref{frac_t2} is as follows.
\begin{quote}
	\begin{verbatim}
	alpha = Variable(1, sign = 'Positive')
	x = Variable(n)
	# specify p and q here
	prob = Problem(Minimize(alpha), [p <= q*alpha])
	\end{verbatim}
\end{quote}

\paragraph{Numerical result.}
Take an example of $p(x) = x^2+1$ and $q(y) = \sqrt{y+0.5}$.
The code for specifying $p$ and $q$ is as the follows.
\begin{quote}
	\begin{verbatim}
p = square(x) + 1
q = sqrt(y+0.5)
	\end{verbatim}
\end{quote}
The global optimal value of the objective function is approximately $1.217$.
With random initial points, DMCP finds the global optimum for problem~\eqref{frac_t1}
 and \eqref{frac_t2}.

\subsection{Linear transceiver design}
\paragraph{Problem description.}
Suppose that a signal $x\in{\mbox{\bf R}}^n$ passes through
a linear pre-coder $A\in{\mbox{\bf R}}^{n\times n}$,
and is transmitted as $Ax$.
Denote the channel matrix as $C\in{\mbox{\bf R}}^{m \times n}$
 and the additive noise as $e\in{\mbox{\bf R}}^m$,
then the received signal is $y = CAx+e$.
The received signal passing through an equalizer $B\in{\mbox{\bf R}}^{n \times m}$
is decoded as $By$.
The problem of determining $A$ and $B$ is called transceiver deign \cite{1254038}.

In this example, we assume that the signal $x$ is binary
and follows IID Bernoulli distribution, and the noise
$e\sim\mathcal{N}(0,\sigma_\mathrm{e}^2 I)$.
Given the channel matrix $C$, the aim is to design $A$ and $B$
such that the mean squared error $\mathbb{E}\|x-By\|_2^2$,
where the mean is taken over $x$ and $e$, is minimized,
and the transmission power is constrained.

An optimization problem is formulated as the following
\[
\begin{array}{ll}
\mbox{minimize} & \frac{1}{2}\|BCA-I\|_F^2 +\sigma_\mathrm{e}^2 \|B\|_F^2 \\
\mbox{subject to} & \|A\|_F \leq p,
\end{array}
\]
where $B$ and $A$ are the variables. The problem is biconvex.

\paragraph{DMCP specification.}
The code can be written as the following.
\begin{quote}
\begin{verbatim}
A = Variable(n,n)
B = Variable(n,m)
sigma_e = Parameter(1)
cost = square(norm(B*C*A-I,'fro'))/2+square(sigma_e)*square(norm(B,'fro'))
prob = Problem(Minmize(cost), [norm(A, 'fro') <= p])
prob.solve(method = 'bcd')
\end{verbatim}
\end{quote}

\begin{figure}
	\centering
	\includegraphics[width=0.6\textwidth]{transceiver.pdf}
	\caption{Linear transceiver design.}
	\label{fig:precoder_eq}
\end{figure}

\paragraph{Numerical result.}
In an experiment, $n = 10$, $m=15$, $p=10$,
and the channel matrix $C$ is a random matrix with IID
normal distribution.
The signal to noise ratio varies, and for each value of $\sigma_\mathrm{e}$,
we try to solve the problem to get a design of $A$ and $B$.
Each design is tested by $1000$ trials with random signal $x$ and noise $e$
generated from the same distributions as the ones in the design.
The method is run without proximal operator and
with initial point generated from the SVD of the channel matrix $C$.
The mean squared error and the averaged bit error rate are in
Figure~\ref{fig:precoder_eq}.

\subsection{Sparse dictionary learning}
\paragraph{Problem description.}
The aim is to find a dictionary $D\in{\mbox{\bf R}}^{m\times n}$
under which the data matrix $X\in{\mbox{\bf R}}^{m\times T}$
can be approximated by sparse coefficients,
i.e., $X \approx DY$
where $Y\in{\mbox{\bf R}}^{n\times T}$ is a sparse matrix \cite{mairal2009online}.

The optimization problem can be formulated as
\[
\begin{array}{ll}
\mbox{minimize} & \frac{1}{2} \|DY-X\|_F^2 + \alpha \|Y\|_1  \\
\mbox{subject to} & \|D\|_F\leq 1,
\end{array}
\]
where the variables are $Y$ and $D$, and $\alpha>0$ is a parameter.
The problem is biconvex.

\paragraph{DMCP specification.}
The code can be written as follows.
\begin{quote}
	\begin{verbatim}
	D = Variable(m,n)
	Y = Variable(n,T)
	alpha = Parameter(sign = 'Positive')
	cost = square(norm(D*Y-X,'fro'))/2+alpha*norm(Y,1)
	prob = Problem(Minimize(cost), [norm(D,'fro') <= 1])
	prob.solve(method = 'bcd')
	\end{verbatim}
\end{quote}
\begin{figure}
	\centering
	\includegraphics[width=0.5\textwidth]{sparse_DL.pdf}
	\caption{Sparse dictionary learning. }
	\label{fig:sparse_DL}
\end{figure}

\paragraph{Numerical result.}
In an experiment, $X$ is a random normal matrix with $m=10$, $n=20$, and $T=20$.
The parameter $\alpha$ is swept from $10^{-5}$ to $1$. For each value of
$\alpha$, the method is called with random initialization,
and the relative approximation error and the cardinality of $Y$ are shown
as a blue dot in Figure \ref{fig:sparse_DL}.

\subsection{Sparse feedback matrix design}
\paragraph{Problem description.}
To design a sparse linear constant output feedback control $u = Ky$
for the system
$$\dot{x} = Ax+Bu, \quad y = Cx,$$
which results in a decay rate $r$ no less than a given threshold
$\theta>0$ in the closed-loop system,
we consider the following optimization problem
\cite{hassibi1999path, boyd1994linear}
\[
\begin{array}{ll}
\mbox{minimize} & \sum_{ij} |K_{ij}|  \\
\mbox{subject to} & P \succeq I,\quad  r \geq \theta \\
& -2r P \succeq (A+BKC)^T P + P(A+BKC),
\end{array}
\]
where $K$, $P$, and $r$ are variables,
and $A$, $B$, $C$, and $\theta$ are given.
The notation $P \succeq I$ means that $P-I$ is semidefinite.
The problem is biconvex with minimal sets of variables to fix $\{P\}$
and $\{K, r\}$.

\paragraph{DMCP specification.}
The code can be the following.
\begin{quote}
	\begin{verbatim}
	P = Variable(n,n)
	K = Variable(m1,m2)
	r = Variable(1)
	cost = norm(K,1)
	constr = [np.eye(n) << P, r >= theta]
	constr += [(A+B*K*C).T*P+P*(A+B*K*C) << -P*r*2]
	prob = Problem(Minimize(cost), constr)
	prob.solve(method = 'bcd')
	\end{verbatim}
\end{quote}
\paragraph{Numerical result.}
An example with $n=m_1=5$, $m_2=4$, $\theta = 0.01$, and
the following data matrices from \cite{hassibi1999path} is tested.
\begin{align*}
A &= \begin{bmatrix}
-2.45 & -0.90 & 1.53 & -1.26 & 1.76 \\
-0.12 & -0.44 & -0.01 & 0.69 & 0.90\\
2.07 & -1.20 & -1.14 & 2.04 & -0.76\\
-0.59 & 0.07 & 2.91 & -4.63 & -1.15\\
-0.74 & -0.23 & -1.19 & -0.06 & -2.52
\end{bmatrix} \\
B & = \begin{bmatrix}
0.81 & -0.79 & 0 & 0 & -0.95\\
-0.34 & -0.50 & 0.06 & 0.22 & 0.92\\
-1.32 & 1.55 & -1.22 & -0.77 & -1.14\\
-2.11 & 0.32 & 0 & -0.83 & 0.59\\
0.31 & -0.19 & -1.09 & 0 & 0
\end{bmatrix} \\
C & = \begin{bmatrix}
0 & 0 & 0.16 & 0 & -1.78 \\
1.23 & -0.38 & 0.75 & -0.38 & 0\\
0.46 & 0 & -0.05 & 0 & 0\\
0 & -0.12 & 0.23 & -0.12 & 1.14\\
\end{bmatrix} \\
\end{align*}
The initial value $P^0$ is an identity matrix, $r^0 = 1$,
and $K^0$ is an matrix with all zeros.
The result is that $r = 0.01$ and
\begin{align*}
K = \begin{bmatrix}
0 &  0.32 &  0 & 0 \\
0 &  -0.46 &  0 &  0 \\
0 & 0 &  0 & 0 \\
0 & 0 & 0 & 0 \\
0 &  0.11 & 0 & 0 \\
\end{bmatrix},
\end{align*}
which is sparse.
The three nonzero entries are in the second column,
so only the second output needs to be fed back.
In the work \cite{hassibi1999path} with decay rate no less than
$0.35$ another sparse feedback matrix is found
with the same cardinality $3$.

\subsection{Bilinear control}
\paragraph{Problem description.}
A discrete time $m$-input bilinear control system is of the following form
\cite{hours2014parametric}
\begin{align*}
x_{t+1} = A x_t + \sum_{i=1}^m u^i_t B^i x_t, \quad t = 1,\ldots,n-1,
\end{align*}
where $u_t = [u^1_t, \ldots, u^m_t]\in\Omega\subseteq{\mbox{\bf R}}^m$
is the input,
and $x_t\in{\mbox{\bf R}}^{d}$ is the system state at time $t$.
In an optimal control problem with fixed initial state,
given system matrices $A, B^i \in{\mbox{\bf R}}^{n\times n}$
and convex objective functions $f$ and $g$,
an optimization problem can be formulated as
\[
\begin{array}{ll}
\mbox{minimize} &f(x) + g(u)\\
\mbox{subject to}
& x_1 = x_\mathrm{ini} \\
& (x_t, u_t) \in \Omega, \quad  t = 1,\ldots,n-1\\
& x_{t+1} = A x_t + \sum_{i=1}^m u^i_t B^i x_t, \quad t = 1,\ldots,n-1,
\end{array}
\]
where $x_t$ and $u_t$ are variables,
and $\Omega$ is a given convex set describing bounds on $u_t$ and $x_t$.
The problem is multi-convex.

As a special case,
a standard model of D.C.-motor is a bilinear system of the following form \cite{Derese1982}
\[
\dot{x} = A_0 x + u A_1 x + bv,
\]
where the derivative is with respect to time, and
\[
x = \begin{bmatrix}
x^1 \\
x^2
\end{bmatrix},\quad
A_0 = \begin{bmatrix}
-1 & 0 \\
0 & -0.1
\end{bmatrix},\quad
A_1 = \begin{bmatrix}
0 & -19 \\
0.1 & 0
\end{bmatrix},\quad
b = \begin{bmatrix}
20 \\
0
\end{bmatrix}.
\]
The correspondence between model and physical variables is that,
$x^1$ is the armature current, $x^2$ is the speed
of rotation, $u$ is the field current, and $v$ is the armature voltage.
For nominal operation,
\[
x^1 = x^2 = u = v = 1.
\]
A control problem is the braking with short-circuited armature ($v_t =0$).
The field current $u$ is controlled
such that the rotation speed decreases to zero as fast as possible,
and that the armature current is not excessively large.
By discretizing over time and taking $10$ samples per second,
the problem can be formulated in the following form
\[
\begin{array}{ll}
\mbox{minimize} & \|x^2\|_2 \\
\mbox{subject to}
&  x_1 = \begin{bmatrix} 1 \\1\end{bmatrix} \\
& \max_{t=1,\ldots, n} |x^1_t| \leq M \\
& (x_{t+1} - x_t)/0.1 = A_0 x_t + u_t A_1 x_t,
\quad t=1,\ldots,n-1,
\end{array}
\]
where $u\in{\mbox{\bf R}}^{n-1}$ and $x_t = [x^1_t, x^2_t]\in{\mbox{\bf R}}^2$
for $t = 1,\ldots,n$ are variables,
and the notation $x^i = [x^i_1, \ldots, x^i_n]$, $i = 1,2$.
The problem is biconvex if we consider
$x = [x_1,\ldots,x_n]\in{\mbox{\bf R}}^{2\times n}$ as one variable.

\paragraph{DMCP specification.}
The code is as the following.
\begin{quote}
\begin{verbatim}
x = Variable(2,n)
u = Variable(n-1)
constr = [x[:,0] == 1, max_entries(abs(x[0,:])) <= M]
for t in range(n-1):
    constr += [x[:,t+1]-x[:,t] == 0.1*(A0*x[:,t]+A1*x[:,t]*u[t])]
prob = Problem(Minimize(norm(x[1,:])), constr)
prob.solve(method = 'bcd')
\end{verbatim}
\end{quote}

\paragraph{Numerical result.}
We take an example with $n=100$ and $M=8$.
The initial value of $x$ is zero
and of $u$ is a vector linearly decreasing from $0.5$ to $0$.
The result is shown in Fig.~\ref{fig:DC_motor},
where the braking is faster than that in a linear control system shown in
\cite{Derese1982}.

\begin{figure}
	\centering
	\includegraphics[width=0.5\textwidth]{DC_motor_1.pdf}
	\caption{D.C.-motor braking.}
	\label{fig:DC_motor}
\end{figure}

\subsection{Resistance estimation}
\paragraph{Problem description.}
A problem in direct current (DC) circuit is to estimate the values of resistors such that
certain constraints on currents and voltages can be satisfied.
The topology of the circuit is given, and several observations
on currents and voltages are known.
A general problem of estimating the resistance to fit the topology and the observations
 can be written as the following
\[
\begin{array}{ll}
\mbox{minimize} & f(u)+g(i) \\
\mbox{subject to} & A(r) i = u \\
& u \in\mathcal{U},\quad  i \in \mathcal{I},\quad r\in\mathcal{R},
\end{array}
\]
where $u\in{\mbox{\bf R}}^n$, $i\in{\mbox{\bf R}}^m$, and $r\in{\mbox{\bf R}}_+^d$ are variables
representing voltages, currents, and resistance, respectively.
The convex functions $f$ and $g$ penalize deviations from the observations.
The first constraint corresponds to the Ohm's law,
where the mapping $A$ is linear and depends on the topology of the circuit.
The sets $\mathcal{U}$, $\mathcal{I}$, and $\mathcal{R}$ are convex,
and they may describe the Kirchhoff's circuit laws.
The problem is multi-convex due to the first constraint.

A simple example is shown in the following circuit diagram.
\begin{center}
\begin{tikzpicture}[scale=1.4]
\draw[color=black]
(-0.2,0)node[]{$u_0$}
(0,0) to [R, l=$a_1$, -*] (0,2)
(1.5,0) to [R, l=$a_2$, *-*] (1.5,2)
(3,0) to [R, l=$a_3$, *-*] (3, 2)
(3,0) to (3.2,0)
(3.75,0) node[]{\large{$\cdots$}}
(4.3,0) to (4.5,0) to [R, l=$a_{n-1}$, *-*] (4.5, 2)
(6,0) to [R, l=$a_n$, *-*] (6, 2)
(7.5,0) to[american current source, l=$I_0$] (7.5,4)

(0,2) to [R, l=$b_1$, *-] (0,4)
(1.5,2) to [R, l=$b_2$, *-*] (1.5,4)
(3,2) to [R, l=$b_3$, *-*] (3,4) to (3.2,4)
(3.75,4) node[]{\large{$\cdots$}}
(4.5,2) to [R, l=$b_{n-1}$, *-*] (4.5,4) to (4.3,4)
(6,2) to [R, l=$b_n$, *-*] (6,4)

(0,2) to [R, l=$c_1$](1.5,2)
(1.5,2) to [R, l=$c_2$,*-*] (3,2) to (3.2,2)
(3.75,2) node[]{\large{$\cdots$}}
(4.3,2) to (4.5,2) to [R, l=$c_{n-1}$,*-*] (6,2)

(-0.2,2.1)node[]{$v_1$}
(1.3,2.1)node[]{$v_2$}
(2.8,2.1)node[]{$v_3$}
(4.2,2.1)node[]{$v_{n-1}$}
(5.8,2.1)node[]{$v_n$}

(0.75,0.2){} node[]{$i_1$}
(2.25,0.2){} node[]{$i_2$}
(5.25,0.2){} node[]{$i_{n-1}$}
(6.75,0.2){} node[]{$i_n$}

(0.75,4.2){} node[]{$j_1$}
(2.25,4.2){} node[]{$j_2$}
(5.25,4.2){} node[]{$j_{n-1}$}
(6.75,4.2){} node[]{$j_n$}

(0,4) to node[currarrow] {} (1.5,4)
(1.5,4) to node[currarrow] {} (3,4)
(4.5,4) to node[currarrow] {} (6,4)
(6,4) to node[currarrow] {} (7.5,4)
(7.5,4) node[ground, rotate = 90]{}

(1.5,0) to node[currarrow,rotate=180] {} (0,0)
(3,0) to node[currarrow,rotate=180] {} (1.5,0)
(6,0) to node[currarrow,rotate=180] {} (4.5,0)
(7.5,0) to node[currarrow,rotate=180] {} (6,0);
\end{tikzpicture}
\end{center}
The known quantities are the
current source $I_0$ and the voltage level $u_0$.
It is observed that $v_k-v_{k+1} \approx \delta$ for $k = 1,\ldots,n-1$,
so the optimization problem is
\[
\begin{array}{ll}
\mbox{minimize} & \sum_{k=1}^{n-1}(v_k-v_{k+1} - \delta)^2\\
\mbox{subject to}
&  x_{1} = y_{1}+z_{1}, \quad x_{n} + z_{n-1} = y_{n}\\
& i_1 = x_1, \quad i_n = -I_0, \quad j_1 = y_1,\quad j_n = -I_0\\
&  x_{k+1} + z_k = y_{k+1}+z_{k+1}, \quad k=1,\ldots,n-2 \\
& x_k a_k = u_0 - v_k, \quad y_k b_k = v_k, \quad k=1,\ldots,n \\
& i_{k+1} = i_k + x_{k+1}, \quad j_{k+1} = j_k + y_{k+1},
\quad z_k c_k = v_k - v_{k+1}, \quad k=1,\ldots,n-1,\\
\end{array}
\]
where $x, y, i, j\in{\mbox{\bf R}}^n, z\in{\mbox{\bf R}}^{n-1}$ are variables for currents,
$a, b\in{\mbox{\bf R}}_+^n, c\in{\mbox{\bf R}}_+^{n-1}$ are the variables for resistance,
and $v\in{\mbox{\bf R}}^n$ is the variable for voltages.
The problem is multi-convex, and the minimal sets to fix are not obvious.

\paragraph{DMCP specification.}
The code can be as the following.
\begin{quote}
\begin{verbatim}
x = Variable(n)
y = Variable(n)
z = Variable(n-1)
i = Variable(n)
j = Variable(n)
a = Variable(n, sign = 'Positive')
b = Variable(n, sign = 'Positive')
c = Variable(n-1, sign = 'Positive')
v = Variable(n)
constr = [x[0] == y[0]+z[0], x[n-1]+z[n-2] == y[n-1]]
constr += [i[0] == x[0], j[0] == y[0], i[n-1] == -I0, j[n-1] == -I0]
cost = 0
for k in range(n-2):
    constr += [x[k+1]+z[k] == y[k+1]+z[k+1]]
for k in range(n):
    constr += [x[k]*a[k] == u0 - v[k], y[k]*b[k] == v[k]]
for k in range(n-1):
    cost += square(v[k]-v[k+1]-delta)
    constr += [z[k]*c[k] == v[k]-v[k+1]]
    constr += [i[k+1]== i[k]+x[k+1], j[k+1] == j[k]+y[k+1]]
prob = Problem(Minimize(cost), constr)
prob.solve(method = 'bcd')
\end{verbatim}
\end{quote}
The \verb|find_minimal_sets| function returns
\begin{quote}
	\begin{verbatim}
[[5, 7, 8], [3, 5, 6], [1, 3, 7], [1, 3, 6],
[1, 7, 8], [5, 6, 8], [1, 6, 8], [3, 5, 7]]
	\end{verbatim}
\end{quote}
where indices $1,3,5,6,7,8$ correspond to variables
$z,a,c,b,y,x$ respectively.

\paragraph{Numerical result.}
We set $n = 10, I_0 = -100, \delta = 1, u_0 = 12$,
and all variables are set with initial value $1$.
The method finds a feasible point with objective value $0$
which solves the problem globally,
and the solution is shown in the following table.

\begin{center}
\begin{tabular}{cccc}
	$a_k$ & $b_k$ & $c_k$ & $v_k-v_{k+1}$ \\
0.120 & 1.123 & 1.800& 1.001 \\
0.224 & 0.987 & 1.776 & 1.002\\
0.329 & 0.885 & 2.249 & 1.003\\
0.431 & 0.785 & 2.904 & 1.004\\
0.535 & 0.687 & 3.653 & 1.004\\
0.641 & 0.591 & 4.066 & 1.004\\
0.743& 0.492 & 3.297 & 1.003\\
0.837 & 0.389 & 2.031 & 1.002\\
0.904 & 0.277 & 1.283 & 1.001\\
0.931& 0.149 &- &-
\end{tabular}
\end{center}

\subsection{Steady state of Markov chain}
\paragraph{Problem description.}
Suppose that $P_1, \ldots, P_n\in{\mbox{\bf R}}^{m\times m}$
are $n$ transition matrices of Markov chains,
then it is known that any convex combination $P=\sum_{i=1}^n \theta_i P_i$ for
$\theta_i \geq 0$ and $\sum_{i=1}^n\theta_i = 1$ is also a transition matrix of
a Markov chain.
Given $P_i$ for $i=1,\ldots,n$ and a convex function $f: {\mbox{\bf R}}^m\to {\mbox{\bf R}}$,
the problem is to find such a convex combination,
so that the Markov chain with respect to $P$ has a steady state vector
$x\in{\mbox{\bf R}}^m$ that achieves the minimum of $f$.

The problem can be formulated as
\[
\begin{array}{ll}
\mbox{minimize} & f(x) \\
\mbox{subject to} & P=\sum_{i=1}^n \theta_i P_i \\
& \theta_i \geq 0, \quad \sum_{i=1}^n\theta_i = 1 \\
& x_j \geq 0, \quad \sum_{j=1}^m x_j = 1 \\
& P^T x =  x,\\
\end{array}
\]
where $x$, $P$, and $\theta$ are variables.
The problem is biconvex.

As an example, $f(x) = \|x-x_0\|_2$, so the goal is
to generate a transition matrix so that
the steady state vector is close to a given distribution $x_0$.

\paragraph{DMCP specification.}
The code is as follows.
\begin{quote}
\begin{verbatim}
cost = norm(x-x0)
constr = [theta >= 0, sum_entries(theta) == 1, x >= 0, sum_entries(x) == 1]
right = 0
for i in range(n):
    right += theta[i]*P0[i]
constr += [P == right, P.T*x == x]
prob = Problem(Minimize(cost), constr)
prob.solve(method = 'bcd')
\end{verbatim}
\end{quote}

\paragraph{Numerical result.}
An example with $n=4$, $m=3$, randomly generated $P_i$, and
\[
x_0 = \begin{bmatrix}
0.25\\
0.3\\
0.45
\end{bmatrix},
\]
is tested.
The initial values are random.
The result gives
\[
x= \begin{bmatrix}
0.25\\
0.3\\
0.45
\end{bmatrix},\quad
\theta = \begin{bmatrix}
0.27\\
0.61\\
0.04\\
0.08
\end{bmatrix},
\]
which achieves the targeted steady state vector.

\subsection{Blind deconvolution}
\paragraph{Problem description.}
Blind deconvolution is an inverse problem commonly encountered in many
practical applications such as image restoration, system identification,
and channel estimation.
The problem is to find two vectors with some priors,
such that their convolution approximates the given data.

Suppose that the data $d\in{\mbox{\bf R}}^{m+n-1}$ is the convolution of
an unknown sparse vector $x_0\in{\mbox{\bf R}}^n$ and an unknown vector $y_0\in{\mbox{\bf R}}^m$.
A problem can be formulated as the following
\[
\begin{array}{ll}
\mbox{minimize} & \|x * y-d\|_2+\alpha \|x\|_1\\
\mbox{subject to} & \|y\|_\infty \leq M,
\end{array}
\]
where $x\in{\mbox{\bf R}}^n$ and $y\in{\mbox{\bf R}}^m$ are variables,
and $\alpha>0$ is a parameter.
The problem is biconvex.
For any $x*y = d$ and a positive scalar $k$
the convolution of $kx$ and $y/k$ is also $d$,
so the constraint $\|y\|_\infty \leq M$ is needed to exclude a trivial solution
$x\approx 0$.

\paragraph{DMCP specification.}
The code can be written as the following.
\begin{quote}
	\begin{verbatim}
	y = Variable(m)
	x = Variable(n)
	cost = norm(conv(y,x)-d,2) + alpha*norm(x,1)
	prob = Problem(Minimize(cost), [norm(y,'inf') <= M])
	prob.solve(method = 'bcd')
	\end{verbatim}
\end{quote}

\paragraph{Numerical result.}
In an example, $m=100$, $n=40$, $M=10$, $\alpha = 0.28$,
and the initial value of every variable is a vector of all ones.
The result is in Figure~\ref{fig:blind_dconv}.

\begin{figure}
	\centering
		\centering
		\includegraphics[height=1.8in]{blind_dconv.pdf}

	\caption{Blind deconvolution.

	\label{fig:blind_dconv}
\end{figure}

\subsection*{Acknowledgments}
This material is based
upon work supported by the National Science Foundation Graduate Research Fellowship
under Grant No. DGE-114747, by the DARPA X-DATA and SIMPLEX programs,
and by the CSC State Scholarship  Fund.
\clearpage


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

\vspace*{1mm}
\noindent
At present the 3-loop heavy flavor corrections form the missing link to analyze the World deep-inelastic scattering
data at next-to-next-to leading order with respect to the determination of the unpolarized parton distribution functions (PDFs)
\cite{PDF},
the strong coupling constant \cite{ALPHA}, and the masses of the charm and bottom quarks \cite{MCMB}. These
distributions and parameters form an essential input for all precision measurements at the LHC, notably the measurement of
the Higgs boson \cite{HIGGS} and $t\bar{t}$ \cite{TTBAR} production cross sections. Furthermore, their precise knowledge
provides new unprecedented tests of the Standard Model and may help in this way to reveal potential deviations pointing to
new physics.

The present project is devoted to calculate the 3-loop heavy flavor Wilson coefficients at large scales $Q^2 \gg m^2$ in analytic
form. In case of the structure function $F_2(x,Q^2)$ this approximation suffices at the $1\%$ level for scales $Q^2/m^2 \raisebox{-0.07cm
10$ \cite{Buza:1995ie}, a region to be selected anyway in case of charm to stay off  higher twist terms \cite{Alekhin:2012ig} at
this precision.

In 2009 a series of Mellin moments has been calculated for the corresponding massive operator matrix elements (OMEs) ranging
from $N=2$ to $N=10...14$, depending on the Wilson coefficient, in Refs.~\cite{Bierenbaum:2009mv} by mapping the operator matrix
elements for the different moments onto massive tadpoles, which were computed using the code {\tt MATAD} \cite{Steinhauser:2000ry}.
These moments serve now
for comparison in the computation of the general-$N$ results. In the time since, four out of five neutral current massive 3-loop
Wilson coefficients have been calculated \cite{Ablinger:2010ty,Behring:2014eya,Ablinger:2014vwa,Ablinger:2014nga,Blumlein:2006mh},
along with seven out of eight massive OMEs
\cite{Ablinger:2010ty,Behring:2014eya,Ablinger:2014vwa,Ablinger:2014nga,Blumlein:2012vq,Ablinger:2014uka,Ablinger:2014lka,
Ablinger:2016eyd,AGG}. For all processes the logarithmic 3-loop contributions are known \cite{Behring:2014eya}. To perform these
computations, several technical,
computer-algebraic and mathematical innovations were necessary, which are described in
Refs.~\cite{Blumlein:1998if,Vermaseren:1998uu,Blumlein:2003gb,Blumlein:2009ta,Blumlein:2009fz,Blumlein:2009cf
Blumlein:2009tj,Blumlein:2010zv,Moch:2001zr,Ablinger:2013cf,Ablinger:2014bra,Ablinger:2011te,
Ablinger:2012qm,Ablinger:2014yaa,Ablinger:2015tua,
Ablinger:2016yjz,SURV}. The different physics results, including also necessary 2-loop calculations, for the neutral
and
charged current reactions  were published in
Refs.~\cite{Bierenbaum:2007qe,Bierenbaum:2008yu,Bierenbaum:2009zt,Blumlein:2014fqa,Blumlein:2016xcy}
and recent surveys were given in \cite{Blumlein:2014zxa,Ablinger:2016eyd}.

In this note we describe recent developments of the project. The paper is organized as follows. The basic formalism and
a series of technical details of the calculation are summarized in Section~\ref{sec:2}. In Section~\ref{sec:3} we review the
status of the calculation of the neutral current structure function $F_2(x,Q^2)$. A survey on the calculation of the
3-loop
corrections to non-singlet charged and neutral current structure functions is given in Section~\ref{sec:4}. Results on the
calculation of the heavy flavor corrections in the full region of $Q^2$ for various  polarized and unpolarized non-singlet structure
functions to $O(\alpha_s^2)$ and associated sum rules are given in Section~\ref{sec:5}. Recent results of 3-loop two-mass corrections
are reported in Section~\ref{sec:6}, and Section~\ref{sec:7} contains the conclusions.

\vspace*{-2mm}
\section{Basic Formalism and Technical Aspects of the Calculation}
\label{sec:2}

\vspace*{1mm}
\noindent
As it has been outlined in Refs.~\cite{Bierenbaum:2009mv,Behring:2014eya} the massive 3-loop Wilson coefficients for
deep-inelastic scattering can be represented in terms of the massless Wilson coefficients \cite{ZN,Vermaseren:2005qc} and
massive operator matrix elements in the asymptotic region $Q^2 \gg m^2$, if the single heavy quark case is considered.
Here $Q^2$ denotes the virtuality of the process and $m$ the heavy quark mass. In the present project the massive
OMEs are newly computed and used in the calculation of the corresponding massive Wilson coefficients. Furthermore, the massive OMEs
describe the transition of a single heavy flavor becoming light in the variable flavor number scheme (VFNS),
cf.~\cite{Buza:1996wv,Behring:2014eya}. Up to 2-loop order the Wilson coefficients are know (semi)analytically
\cite{Laenen:1992zk}\footnote{For a precise implementation in Mellin space, see \cite{Alekhin:2003ev}.}. In the asymptotic region
analytic
expressions have been derived in Refs.~\cite{Buza:1995ie,Bierenbaum:2007qe}.

The whole 3-loop project constitutes of a calculation of 2864 Feynman diagrams mapping to some dozens of $10^4$ scalar integrals.
They are generated using {\tt QGRAF} \cite{QGRAF}, calculating their Dirac structure using {\tt FORM} \cite{FORM} and the color
structure using {\tt Color} \cite{vanRitbergen:1998pn}.
Except for simpler topologies, we use the integration by part (IBP)-reduction \cite{IBP} implemented in the package {\tt Reduze\;2}
\cite{REDUZE}\footnote{The package uses the codes {\tt GiNaC} \cite{Bauer:2000cp} and {\tt Fermat} \cite{FERMAT}.} to map
the major part of the problem to the calculation of 687 master integrals, out of which 116 remain to be calculated. 571 master integrals
can be solved by sum-structures using difference-field theory
\cite{Karr:81,Schneider:01,Schneider:05a,Schneider:07d,Schneider:10b,Schneider:10c,Schneider:15a,Schneider:08c,DFTheory} implemented in
the
package {\tt Sigma} \cite{SIG1,SIG2}, or other techniques, like  the use of hypergeometric functions \cite{HYP}, Mellin-Barnes
representations \cite{MB}, the method of hyperlogarithms \cite{HYPLOG,Ablinger:2014yaa}, the solution of differential
equations \cite{DEQ}, and the Almkvist-Zeilberger algorithm \cite{AZ},
implemented within the package {\tt MultiIntegrate} \cite{Ablinger:PhDThesis}.

They can be expressed by iterative sum-structures in $N$-space or iterative integrals in $x$-space, i.e. in terms of harmonic
polylogarithms
\cite{Remiddi:1999ew}, generalized harmonic polylogarithms of the Kummer-type \cite{Moch:2001zr,Ablinger:2013cf}, cyclotomic harmonic
polylogarithms \cite{Ablinger:2011te}, or root-valued
iterated integrals \cite{Ablinger:2014bra}, which all correspond to difference or differential equations that can be factorized to
first order.
This is not the case for all the remaining 116 master integrals, which depend on structures obeying 2nd order equations. Here we
expect complete elliptic integrals of rational argument  and related functions to emerge, see \cite{ELLIPT}\footnote{Similar structures
were observed also in Refs.~\cite{ERLT,Adams:2016xah,Adams:2015ydq}.}, over which first order structures are iterated again. We
currently work on the solution of these systems.

In all these methods the solution of recurrences and the various properties of special functions emerging in this context play a central
role and have to be used algorithmically. This is made possible by the packages {\tt Sigma} \cite{SIG1,SIG2}, {\tt EvaluateMultiSums} and
{\tt SumProduction}, \cite{EMSSP}, {\tt RhoSum} \cite{RHOSUM}, decoupling formalisms \cite{DECOUP}, {\tt
HarmonicSums} \cite{Ablinger:PhDThesis,Harmonicsums,Ablinger:2013cf,Ablinger:2011te,Ablinger:2014bra}
and {\tt MultiIntegr- ate} \cite{Ablinger:PhDThesis}.

\vspace*{-2mm}
\section{Status of the 3-Loop Neutral Current Corrections}
\label{sec:3}

\vspace*{1mm}
\noindent
The heavy flavor contributions to the structure functions $F_{2,L}(x,Q^2)$ are given by
\begin{eqnarray}
F_{2,L}^{\rm heavy}(x,Q^2) =F_{2,L}(x,Q^2) - F_{2,L}^{\rm massless}(x,Q^2)~,
\end{eqnarray}
for a single heavy quark $Q$ and $N_F$ massless quarks
\begin{eqnarray}
\label{eqF2}
F_{(2,L)}^{\rm heavy}(x,N_F\!\!\!&+&\!\!\!1,Q^2,m^2) =\nonumber\\
       &&\hspace*{-1cm} \sum_{k=1}^{N_F}e_k^2\Biggl\{
                   L_{q,(2,L)}^{\sf NS}\left(x,N_F+1,\frac{Q^2}{m^2}
                                                ,\frac{m^2}{\mu^2}\right)
                \otimes
                   \Bigl[f_k(x,\mu^2,N_F)+f_{\overline{k}}(x,\mu^2,N_F)\Bigr]
\nonumber\\ &&\hspace{4mm}
               +\frac{1}{N_F}L_{q,(2,L)}^{\sf PS}\left(x,N_F+1,\frac{Q^2}{m^2}
                                                ,\frac{m^2}{\mu^2}\right)
                \otimes
                   \Sigma(x,\mu^2,N_F)
\nonumber\\ &&\hspace{4mm}
               +\frac{1}{N_F}L_{g,(2,L)}^{\sf S}\left(x,N_F+1,\frac{Q^2}{m^2}
                                                 ,\frac{m^2}{\mu^2}\right)
                \otimes
                   G(x,\mu^2,N_F)
                             \Biggr\}
\nonumber\\
       &+&  e_Q^2\Biggl[
                   H_{q,(2,L)}^{\sf PS}\left(x,N_F+1,\frac{Q^2}{m^2}
                                        ,\frac{m^2}{\mu^2}\right)
                \otimes
                   \Sigma(x,\mu^2,N_F)
\nonumber\\ && \hspace{7mm}
                  +H_{g,(2,L)}^{\sf S}\left(x,N_F+1,\frac{Q^2}{m^2}
                                           ,\frac{m^2}{\mu^2}\right)
                \otimes
                   G(x,\mu^2,N_F)
                                  \Biggr]~.
\end{eqnarray}
Here $e_k$ and $e_Q$ denote the light and heavy quark charges, $f_k$ and $G$ are the quarkonic and gluon
parton densities and $\Sigma = \sum_{k=1}^{N_F} \left(f_k + f_{\bar{k}}\right)$ the singlet density,
and $\mu$ denotes the factorization
scale. The five heavy flavor Wilson coefficients are given by $L_{q,}^{\sf NS}, L_{q}^{\sf PS}, L_{g}^{\sf S}$ and
$H_q^{\sf PS}, H_g^{\sf S}$, cf.~Refs.~\cite{Ablinger:2010ty,Behring:2014eya,Ablinger:2014vwa,Ablinger:2014nga}.

In Figure~\ref{fig:1} the present status of the 3-loop results on $F_2^{\rm charm}(x,Q^2)$ is illustrated for $Q^2 = 100~{\rm GeV}^2$
referring to the parton densities \cite{Alekhin:2013nda} and $m_c^{\rm pole} = 1.59~{\rm GeV}$. We show the different contributions
in $O(\alpha_s), O(\alpha_s^2)$ and $O(\alpha_s^3)$ using the asymptotic representation,
highlighting their behaviour at large values of $x$ in the inset.

\noindent
\begin{figure}\centering
     \includegraphics[width=0.7\textwidth]{F2-defense-Lq2PS-as3.pdf}
\caption[]{The different known contributions to $F_{2}^{\rm charm}(x,Q^2)$ by the five heavy flavor Wilson coefficients
from $O(\alpha_s)$ (dotted lines),  $O(\alpha_s^2)$ (dashed lines), $O(\alpha_s^2)$ (full lines),
for $m_c^{\rm pole} = 1.59~{\rm GeV}$ and the PDFs \cite{Alekhin:2013nda}. The contributions due to $H_g^{\sf S}$ have been scaled down by
a factor of 20 for better visibility.
\label{fig:1}
}
\end{figure}

\vspace*{-3mm}
The intended accuracy is $O(1\%)$ and requires all the contributions mentioned above, at least in some part of the
kinematic range.
The largest contribution is due to $H_g^{\sf S}$ driven by the gluon distribution starting at $O(\alpha_s)$. Here the calculation of the
3-loop term is underway, while all other contributions are known already. The next largest term in the small $x$ region is
$H_q^{\sf
PS}$, yielding negative corrections there. It is followed by $L_g^{\sf S}$ with smaller positive corrections at 3-loop order, which are
larger than those in 2-loop order by the same Wilson coefficient. The fourth largest contribution is due to $L_q^{\sf PS}$ and
the smallest contributions are due to $L_q^{\sf NS}$, while it is largest in the valence region. The full knowledge of the 3-loop heavy
flavor corrections to $F_2(x,Q^2)$ will
improve on presently remaining theory errors both on $\alpha_s(M_Z^2)$ and $m_c$ \cite{MCMB} due to a yet approximate treatment
\cite{Kawamura:2012cr} of these corrections, based on the previous work
\cite{Bierenbaum:2009mv,Bierenbaum:2007qe,Bierenbaum:2008yu,Ablinger:2010ty,Bierenbaum:2009zt}.

\vspace*{-2mm}
\section{3-Loop Non-Singlet Polarized and Charged Current Corrections}
\label{sec:4}

\vspace*{1mm}
\noindent
The flavor non-singlet OME \cite{Ablinger:2014vwa} appears in various unpolarized and polarized neutral
and charged current structure functions, combined with the different massless Wilson coefficients and helicity
dependent anomalous dimensions \cite{Vermaseren:2005qc,Moch:2008fj,Davies:2016ruz,Moch:2014sna}. This allows to
calculate the polarized 3-loop massive non-singlet Wilson coefficient for twist-2 part of the structure
functions $g_{1,2}(x,Q^2)$ and the charged current structure functions $xF_3^{W^+ - W^-}(x,Q^2)$ and
$F_{1,2}^{W^+ - W^-}(x,Q^2)$, cf.~Refs.~\cite{Behring:2015zaa,Behring:2015roa,F1F2}; for the lower order corrections
see \cite{Blumlein:2014fqa,CC}. Here, $g_2$ is obtained
by the Wandzura-Wilczek relation \cite{Wandzura:1977qf,Jackson:1989ph,Blumlein:1996vs,Blumlein:1998nv}.

In Figure~\ref{fig:2} we show the relative 3-loop corrections to $g_{1}^{\rm NS}(x,Q^2)$ and $xF_3^{W^+ +
W^-}(x,Q^2)$ for typical scales $Q^2$ using the polarized PDFs \cite{Blumlein:2010rn} and the
unpolarized PDFs \cite{Alekhin:2013nda}, normalizing to the massless contributions.
\noindent
\begin{figure}\centering
\includegraphics[width=0.49\textwidth]{g1-ratio-as3.pdf}
\includegraphics[width=0.49\textwidth]{F3-ratio-as3.pdf}
\caption[]{The relative 3-loop charm quark corrections to the polarized neutral current structure function
$g_1^{\rm NS}(x,Q^2)$ using the PDFs \cite{Blumlein:2010rn} and the unpolarized charged current
structure function $xF_3^{W^+ + W^-}(x,Q^2)$ using the
PDFs \cite{Alekhin:2013nda} for $m_c^{\rm pole} = 1.59~{\rm GeV}$ as functions of $x$ and $Q^2$;
from \cite{Behring:2015zaa} and \cite{Behring:2015roa}.
\label{fig:2}}
\end{figure}

\vspace*{-3mm}
The 3-loop heavy flavor corrections to $g_1^{\rm NS}$ vary between $+2\%$ and $-5\%$ and those
for the charged current structure function $xF_3^{W^+ + W^-}$ between $+3\%$ and $-3\%$, compared to
the light flavor contributions. These corrections cannot be resolved in present measurements but
will play a role in high luminosity measurements at planned future colliders \cite{FUTURE}.

The asymptotic 3-loop heavy flavor corrections to the charged current structure functions $F_{1(2)}^{W^+ - W^-}$
is shown in Figure~\ref{fig:3}. The relative contribution ranges from $-1 (-2)\%$ to $-8\%$, with the largest corrections
around $x \sim 0.03$. The higher order terms yield positive corrections.

\vspace*{-2mm}
\section[]{\boldmath Power Corrections in the Non-Singlet Case at
\boldmath {\boldmath $\mathbf{O(\alpha_s^2)}$}}
\label{sec:5}

\vspace*{1mm}
\noindent
In the foregoing sections we have presented results calculating the massive Wilson coefficients in the asymptotic region
$Q^2 \gg m^2$. It is, however, also interesting to study the Wilson coefficients
at lower scales $Q^2$ for data in this
region and to estimate the scales at which the asymptotic representation is valid. We investigate the non-singlet
massive Wilson coefficients to $O(\alpha_s^2)$ for the inclusive structure functions \cite{Blumlein:2016xcy}.

\noindent
\begin{figure}\centering
\includegraphics[width=0.49\textwidth]{F1ratios100asy.pdf}
\includegraphics[width=0.49\textwidth]{F2ratios100asy.pdf}
\caption[]{The relative 3-loop charm quark corrections to the charged current structure functions
$F_{1(2)}^{W^+ - W^-}(x,Q^2)$ at $Q^2 = 100~{\rm GeV}^2$ for $m^{\rm pole} = 1.59~{\rm GeV}$ using the PDFs
of Ref.~\cite{Alekhin:2013nda}; form \cite{F1F2}.
\label{fig:3}}
\end{figure}

\vspace*{-3mm}
\noindent
\begin{figure}\centering
\includegraphics[width=0.49\textwidth]{pBj2.pdf}
\includegraphics[width=0.49\textwidth]{CuBJ2.pdf}
\caption[]{Left panel:~The charm and bottom quark contributions to the polarized Bjorken sum rule as a function of
$\xi_c = Q^2/m_c^2$, for $m_c^{\rm pole} = 1.59~{\rm GeV}$ and $m_b^{\rm pole} = 4.78~{\rm GeV}$ \cite{Agashe:2014kda}.
The function $C_{\rm pBj}$ describes
the flavor excitation, with $C_{\rm pBj} = 1$ one heavy quark being effectively massless. Right panel:~The 2-loop
charm quark contribution to the unpolarized Bjorken sum rule; from \cite{Blumlein:2016xcy}.
\label{fig:4}}
\end{figure}

\vspace*{-3mm}
In the tagged flavor case a similar analysis was performed in Ref.~\cite{Blumlein:1998sh}\footnote{For the
$O(\alpha_s^2)$ results on the tagged contributions to the structure functions see
\cite{Buza:1995ie,Bierenbaum:2007qe}.}. It turns out that the tagged
flavor case does not lead to a stable description at large scales $Q^2$, containing logarithmic terms, which are absent
in the inclusive case due to heavy flavor loop effects in diagrams with massless final states. In the inclusive case
the decoupling of heavy flavors is uniquely described. Going to very low scales $Q^2$ for deep-inelastic scattering
logarithmic terms remain. However, the deep-inelastic description is limited to virtualities $Q^2$ of a certain size,
say $Q^2 \raisebox{-0.07cm    Q_0^2 = 5~{\rm GeV}^2$ and cannot be applied below.

We show the effect of these contributions on the structure functions $F_{1(2)}$ in Figure~\ref{fig:3} (dash-dotted lines)
replacing the contributions up to 2-loops by the complete results. At $Q^2 = 100~{\rm GeV}^2$ the corrections do widely agree
with the asymptotic result but a small difference in the large $x$ range still remains.

One also obtains interesting results on different sum rules for the structure functions, such as
the Adler- \cite{Adler:1965ty}, unpolarized Bjorken- \cite{Bjorken:1967px}, polarized Bjorken- \cite{Bjorken:1969mm},
and Gross-Lewellyn Smith sum rule \cite{Gross:1969jf}. While due to fermion number conservation the contribution to these sum
rules is free of mass effects in the asymptotic range, the last three sum rules receive heavy quark mass effects. We have shown
that the Adler sum rule does not receive any correction, also no target mass corrections.

Figure~\ref{fig:4} demonstrates the transition of the massive quark contributions to massless ones
for the massive 2-loop corrections of the polarized and unpolarized Bjorken sum-rule as a function of $Q^2$.
One should note that the respective corrections enter the $N_F$-terms of the sum rule only, cf.~\cite{Blumlein:2016xcy}.
The transition of this term to 1 only proceeds slowly both for charm and bottom quarks as a function of $Q^2$. At lower
scales $Q^2$ the corresponding contributions are even negative in case of the polarized Bjorken sum rule due to virtual
corrections, which is of importance in experimental analyses.

\vspace*{-2mm}
\section{3-Loop Two-Mass Corrections in the Asymptotic Region}
\label{sec:6}

\vspace*{1mm}
\noindent
Starting at 3-loop order the massive Wilson coefficients contain Feynman diagrams in which two internal
massive fermion lines are present leading to contributions of charm and bottom in individual terms, which
cannot be separated. Due to this the variable flavor number scheme presented  in \cite{Buza:1996wv,Behring:2014eya}
is no longer applicable in the strict sense, but needs to be generalized. As $m_b^2/m_c^2 \sim 10$, one possibility
for large scales $Q^2 \gg m_c^2, m_b^2$ consists in decoupling both heavy quarks together. The corresponding
variable flavor number scheme has been worked out in Ref.~\cite{TWOM2} and requires to obtain as well the
two-mass contributions to the different OMEs at 3-loop order. Their moments $N= 2,4,6$ have already been calculated
\cite{TWOM2,TWOM1} projecting them to massive tadpoles and using the package {\tt Q2e} \cite{Q2E}. In the two-mass
case the renormalization of the OMEs, given in the single mass case in \cite{Bierenbaum:2009mv}, has to be extended.
Furthermore,
the 2-mass OMEs $A_{qq,Q}^{\rm NS},~A_{qq,Q}^{\rm NS,TR}$ and $A_{gq,Q}$ have been calculated for general mass
assignment analytically and the calculation of $A_{Qq}^{\rm PS}$ is underway. In case of the OME $A_{gg,Q}$
all scalar topologies have been calculated. Here new classes of iterative integrals do emerge, which can be handled
with the package {\tt HarmonicSums.m} by now. While an expansion in the mass ratio $m_c^2/m_b^2$ is possible
for fixed moments, it turns out not to be possible for a general Mellin variable $N$, thus requiring
the complete calculation.

\vspace*{-2mm}
\section{Conclusions}
\label{sec:7}

\vspace*{1mm}
\noindent
After having calculated a set of 3-loop moments for all massive OMEs contributing to the massive 3-loop
Wilson coefficients in the asymptotic region $Q^2 \gg m^2$ in 2009 \cite{Bierenbaum:2009mv}, progress
has been made in the calculation of the  general $N$ results of these Wilson coefficients. All logarithmic
contributions are known \cite{Behring:2014eya} and the Wilson coefficients $L_{q,}^{\sf NS}, L_{q}^{\sf PS},
L_{g}^{\sf S}$ and $H_q^{\sf PS}$ have been calculated along with the OMEs appearing as the transition matrix
elements in the variable flavor number scheme at 3-loop order,
cf.~\cite{Ablinger:2010ty,Behring:2014eya,Ablinger:2014vwa,Ablinger:2014nga}.
In case of the charged current structure functions and the non-singlet  contribution to the polarized structure functions
$g_{1(2)}(x,Q^2)$ a series of 3-loop Wilson coefficients have been calculated \cite{Behring:2015zaa,Behring:2015roa,F1F2},
completing the 2-loop charged current programme \cite{CC,Blumlein:2014fqa} before. Recently, the non-singlet $O(\alpha_s^2)$
corrections have been completed for the most important inclusive polarized and unpolarized neutral and charged current structure
functions and the sum rules associated to them in the whole range $Q^2$ relevant for deep-inelastic scattering \cite{Blumlein:2016xcy}.
It has been shown that the process of a heavy
quark becoming effectively massless at high scales proceeds slowly, which should be considered in the matching in the
different schemes being used in data analysis currently. Clearly, the matching at $m_Q^2 = Q^2$ is unfortunate since a heavy quark's
velocity is not ultra-relativistic there.

At 3-loop order two-heavy-mass Feynman diagrams contribute to the OMEs. This requires a different renormalization compared to that of
the single-heavy-mass case and leads to a change of the associated variable flavor number scheme since these terms
are neither charm nor bottom contributions. One may, however, design a variable flavor number scheme in decoupling both contributions
together. Beyond a series of moments also the general $N$ contributions have been computed in some cases \cite{TWOM2}.

We have devised an algorithm to calculate massive Feynman diagrams containing local operator insertions, mapping the
differential equations obtained from the IBP-relations into systems of difference equations, which are solved automatically
\cite{Ablinger:2015tua,Ablinger:2016yjz}. This algorithm factorizes the corresponding problem to first order structures as far as
possible and into
potential remaining terms, which cannot be factorized neither in $N$- nor $x$-space. If the latter terms are not present,
we receive an iterative sum or integral solution of the corresponding physical problem starting from differential equations in whatever
basis and obtaining the emerging letter-representation of the contributing alphabet in an automatic manner together
with a proof
certificate. In this way we automatically detect non-factorizing contributions uniquely. They need a further separate
treatment to obtain the solution.

In the present project various new mathematical structures have been found for Feynman-integrals in general.
While at the time of 1998 the harmonic sums yielded a sufficient systematic representation for massless 2-loop
problems \cite{Blumlein:1998if,Vermaseren:1998uu}, the present massive 3-loop problems require Kummer-iterated
integrals \cite{Moch:2001zr,Ablinger:2013cf}, cyclotomic iterated integrals \cite{Ablinger:2011te}, root-valued iterated
integrals \cite{Ablinger:2014bra}, and further generalizations to elliptic and $_2F_1$-valued non-iterative letters in
otherwise iterative integrals \cite{ELLIPT}. This list of structures will be growing further addressing processes at
even higher loops and with more scales. The packages {\tt Sigma, EvaluateMultiSums, SumProduction, RhoSum}, {\tt
HarmonicSums} and {\tt MultiIntegrate}
\cite{SIG1,SIG2,EMSSP,RHOSUM,DECOUP,Ablinger:PhDThesis,Harmonicsums,Ablinger:2013cf,Ablinger:2011te,Ablinger:2014bra}
were significantly extended or newly created in the project and are already in partly use in a series of other projects.

Presently we work on the completion of $A_{Qg}^{(3)}$, in which new mathematical structures appear, the solution of which has to be
automated.

\vspace*{-2mm}