# Continuants and convergence of certain continued fractions

Daniel Duverney  110, rue du chevalier français, 59000 Lille, France
<EMAIL_ADDRESS>and  Iekata Shiokawa  13-43, Fujizuka-cho, Hodogaya-
ku, Yokohama 240-0031, Japan<EMAIL_ADDRESS>
(Date: October 18, 2022)

###### Abstract.

We give a concise introduction to the theory of continuants and show how
Perron used them in his proof of Tietze theorem on the convergence of infinite
semi-regular continued fractions, as well as for the study of the convergence
of purely periodic continued fractions.

###### Key words and phrases:

Continuant, determinant, Tietze Theorem, Stolz Theorem, Convergence, Semi-
regular continued fractions, Purely periodic continued fractions.

###### 1991 Mathematics Subject Classification:

11A55, 15A69

##  1\. Introduction

The purpose of this note is to give a short introduction to continuants and
show how they can be used in order to study the convergence of continued
fractions in some cases. Continuants and their relations to continued
fractions were first considered by Spottiswoode in 1856 [12]. A few years
later, Nachreiner [8] and Muir [7] studied them independently. Muir called
them continuants for the first time.

First we recall some basic facts about continued fractions [5]. Let $a_{1},$
$a_{2},\ldots,$ $b_{0},$ $b_{1},$ $b_{2},\ldots$ be infinite sequences of
indeterminates. We define

(1.1) $A_{-1}=1,\quad A_{0}=b_{0},\quad B_{-1}=0,\quad B_{0}=1,$

and for $n\geq 1$

(1.2) $\displaystyle A_{n}$ $\displaystyle=b_{n}A_{n-1}+a_{n}A_{n-2},$ (1.3)
$\displaystyle B_{n}$ $\displaystyle=b_{n}B_{n-1}+a_{n}B_{n-2}.$

It is well known that

(1.4)
$\alpha_{n}:=b_{0}+\frac{a_{1}}{b_{1}}\genfrac{}{}{0.0pt}{}{{}}{+}\frac{a_{2}}{b_{2}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots+}\frac{a_{n}}{b_{n}}=\frac{A_{n}}{B_{n}}\qquad\left(n\geq
1\right),$

and an easy induction using (1.2) and (1.3) shows that

(1.5) $A_{n}B_{n-1}-A_{n-1}B_{n}=\left(-1\right)^{n-1}a_{1}a_{2}\cdots
a_{n}\qquad\left(n\geq 1\right).$

This yields immediately

(1.6)
$\frac{A_{n}}{B_{n}}-\frac{A_{n-1}}{B_{n-1}}=\left(-1\right)^{n-1}\frac{a_{1}a_{2}\cdots
a_{n}}{B_{n-1}B_{n}}\qquad\left(n\geq 1\right).$

Continuants allow to generalize (1.2) and (1.3) by expressing $A_{n+k}$ and
$B_{n+k}$ in terms of $A_{n-1},$ $A_{n-2},$ $B_{n-1},$ $B_{n-2}$ for all
$k\geq 0$. This leads consequently to a generalization of (1.5) and (1.6), as
we will see in Section 2, which consists of a short introduction to
continuants.

In Section 3, we will show how Perron used continuants in [10, Chapter 5] for
proving Theorem 1 below, known as Tietze theorem. By Theorem 1, any infinite
semi-regular continued fraction ([15],[10],[2]) is convergent. An alternative
proof of Theorem 1 (not using continuants) can be found in [3].

######  Theorem 1.

Assume that the infinite continued fraction

(1.7)
$\alpha:=b_{0}+\frac{a_{1}}{b_{1}}\genfrac{}{}{0.0pt}{}{{}}{+}\frac{a_{2}}{b_{2}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots+}\frac{a_{n}}{b_{n}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots}$

satisfies the following conditions:

(1.8) $a_{n}\in\left\\{-1,1\right\\},\qquad
b_{n}\in\left[1,+\infty\right[,\qquad b_{n}+a_{n+1}\geq 1\qquad\left(n\geq
1\right).$

Then $\alpha$ is convergent.

In Section 4, we will use continuants for proving the following result on the
convergence of purely periodic continued fractions.

######  Theorem 2.

Let $\left(a_{n}\right)_{n\geq 1}$ and $(b_{n})_{n\geq 0}$ be non-zero complex
numbers. Assume that the infinite continued fraction

$\alpha:=b_{0}+\frac{a_{1}}{b_{1}}\genfrac{}{}{0.0pt}{}{{}}{+}\frac{a_{2}}{b_{2}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots+}\frac{a_{n}}{b_{n}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots}$

is purely periodic with period $p\geq 1,$ that is

$a_{n+p}=a_{n}\quad\left(n\geq 1\right),\qquad b_{n+p}=b_{n}\quad\left(n\geq
0\right).$

Let $\lambda_{1}$ and $\lambda_{2}$ with
$\left|\lambda_{1}\right|\geq\left|\lambda_{2}\right|$ be the eigenvalues of
the matrix

$M:=\left(\begin{array}[c]{cc}A_{p-1}&a_{p}A_{p-2}\\\
B_{p-1}&a_{p}B_{p-2}\end{array}\right)=\left(\begin{array}[c]{cc}A_{p-1}&A_{p}-b_{0}A_{p-1}\\\
B_{p-1}&B_{p}-b_{0}B_{p-1}\end{array}\right).$

Then $\alpha$ is convergent if and only if $B_{p-1}\neq 0$ and one of the
following conditions holds:

(C1) $\lambda_{1}=\lambda_{2}.$

(C2) $\left|\lambda_{1}\right|>\left|\lambda_{2}\right|$ and
$A_{q}-x_{2}B_{q}\neq 0$ for $0\leq q\leq p-2,$ where $x_{2}$ is defined by

(1.9) $\lambda_{2}=x_{2}B_{p-1}+a_{p}B_{p-2}.$

Moreover, in case of convergence, $\alpha=x_{1},$ where $x_{1}$ is defined by

(1.10) $\lambda_{1}=x_{1}B_{p-1}+a_{p}B_{p-2}.$

Theorem 2 was first proved by Stolz [13]. The basic idea of the proof we give
in Section 4 is due to Perron [11, Page 83], although Perron doesn’t use
matrix calculations. It consist in computing the values of $A_{n}$ and $B_{n}$
for all $n\geq 0$ in function of $\lambda_{1}$ and $\lambda_{2}.$ An
alternative treatment, using the properties of linear fractional
transformations, can be found in [5, Section 3.2].

Finally, in Section 5, we will deduce from Theorem 2 a short proof of Theorem
3 below, known as Galois generalized theorem [5, Theorem 3.4]. Recall that the
original Galois theorem applies to regular continued fractions ([4],[10, Satz
3.6], [1, Exercise 4.6]). Let

(1.11)
$\alpha:=b_{0}+\frac{a_{1}}{b_{1}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots+}\frac{a_{p-1}}{b_{p-1}}\genfrac{}{}{0.0pt}{}{{}}{+}\frac{a_{p}}{b_{0}}\genfrac{}{}{0.0pt}{}{{}}{+}\frac{a_{1}}{b_{1}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots+}\frac{a_{p-1}}{b_{p-1}}\genfrac{}{}{0.0pt}{}{{}}{+}\frac{a_{p}}{b_{0}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots}$

be a purely periodic continued fraction, and let

(1.12)
$\alpha^{\prime}:=b_{0}+\frac{a_{p}}{b_{p-1}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots+}\frac{a_{2}}{b_{1}}\genfrac{}{}{0.0pt}{}{{}}{+}\frac{a_{1}}{b_{0}}\genfrac{}{}{0.0pt}{}{{}}{+}\frac{a_{p}}{b_{p-1}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots+}\frac{a_{2}}{b_{1}}\genfrac{}{}{0.0pt}{}{{}}{+}\frac{a_{1}}{b_{0}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots},$

its reverse continued fraction. Note that $\alpha^{\prime}$ is purely
periodic. Define $a_{n}^{\prime}$ and $b_{n}^{\prime}$ $(n\geq 1)$ by

$\alpha^{\prime}=b_{0}^{\prime}+\frac{a_{1}^{\prime}}{b_{1}^{\prime}}\genfrac{}{}{0.0pt}{}{{}}{+}\frac{a_{2}^{\prime}}{b_{2}^{\prime}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots+}\frac{a_{n}^{\prime}}{b_{n}^{\prime}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots}$

and $A_{n}^{\prime},$ $B_{n}^{\prime}$ $(n\geq-1)$ by $A_{-1}^{\prime}=1,$
$A_{0}^{\prime}=b_{0}^{\prime}=b_{0},$ $B_{-1}^{\prime}=0,$
$A_{-1}^{\prime}=1$ and

$A_{n}^{\prime}=b_{n}^{\prime}A_{n-1}^{\prime}+a_{n}^{\prime}A_{n-2}^{\prime},\quad
B_{n}^{\prime}=b_{n}^{\prime}B_{n-1}^{\prime}+a_{n}^{\prime}B_{n-2}^{\prime}\quad\left(n\geq
1\right).$

######  Theorem 3.

Assume that the purely periodic continued fraction (1.11) is convergent, and
let $\lambda_{1},$ $\lambda_{2},$ $x_{1},$ and $x_{2}$ be as in Theorem 2.
Then $\alpha=x_{1}$ and its reverse continued fraction (1.12) converges to
$b_{0}-x_{2}$, except if $\left|\lambda_{1}\right|>\left|\lambda_{2}\right|$
and there exists $q\in\left\\{0,\ldots,p-2\right\\}$ such that

(1.13) $A_{q}^{\prime}-\left(b_{0}-x_{1}\right)B_{q}^{\prime}=0$

in which case it is divergent.

######  Corollary 1.

Assume that the $a_{n}$ $\left(n\geq 1\right)$ and $b_{n}$ $(n\geq 0)$ are
non-zero integers and that the continued fraction (1.11) converges to an
irrational $\alpha$111The example
$\alpha:=2-\frac{1}{2}\genfrac{}{}{0.0pt}{}{{}}{-}\frac{1}{2}\genfrac{}{}{0.0pt}{}{{}}{-\cdots-}\frac{1}{2}\genfrac{}{}{0.0pt}{}{{}}{-\cdots}=1$
shows that $\alpha$ can be rational.. Then $\alpha$ is quadratic and

$-\alpha^{\ast}=\alpha^{\prime}-b_{0}=\frac{a_{p}}{b_{p-1}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots+}\frac{a_{2}}{b_{1}}\genfrac{}{}{0.0pt}{}{{}}{+}\frac{a_{1}}{b_{0}}\genfrac{}{}{0.0pt}{}{{}}{+}\frac{a_{p}}{b_{p-1}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots+}\frac{a_{2}}{b_{1}}\genfrac{}{}{0.0pt}{}{{}}{+}\frac{a_{1}}{b_{0}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots},$

where $\alpha^{\ast}$ is the conjugate of $\alpha.$

Corollary 1 will also be proved in Section 5. It applies, for example, to any
purely periodic semi-regular continued fraction $\alpha.$ In the case of a
regular continued fraction, i.e. $a_{n}=1$ for all $n\geq 1,$

$-\frac{1}{\alpha^{\ast}}=b_{p-1}+\frac{1}{b_{p-2}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots+}\frac{1}{b_{1}}\genfrac{}{}{0.0pt}{}{{}}{+}\frac{1}{b_{0}}\genfrac{}{}{0.0pt}{}{{}}{+}\frac{1}{b_{p-1}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots+}\frac{1}{b_{1}}\genfrac{}{}{0.0pt}{}{{}}{+}\frac{1}{b_{0}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots}.$

This is Galois Theorem. Similarly, in the case of a negative continued
fraction, i.e. $a_{n}=-1$ for all $n\geq 1$,

$\frac{1}{\alpha^{\ast}}=b_{p-1}-\frac{1}{b_{p-2}}\genfrac{}{}{0.0pt}{}{{}}{-\cdots-}\frac{1}{b_{1}}\genfrac{}{}{0.0pt}{}{{}}{-}\frac{1}{b_{0}}\genfrac{}{}{0.0pt}{}{{}}{-}\frac{1}{b_{p-1}}\genfrac{}{}{0.0pt}{}{{}}{-\cdots-}\frac{1}{b_{1}}\genfrac{}{}{0.0pt}{}{{}}{-}\frac{1}{b_{0}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots},$

which has been proved by Möbius ([6],[16, Satz 19]).

##  2\. Continuants

We call continuant [7] any determinant of the form

(2.1)
$K\binom{a_{1},\ldots,a_{n}}{b_{0},\ldots,b_{n}}:=\left|\begin{array}[c]{ccccc}b_{0}&-1&0&\cdots&0\\\
a_{1}&b_{1}&-1&\ddots&\vdots\\\ 0&a_{2}&b_{2}&\ddots&0\\\
\vdots&&\ddots&\ddots&-1\\\
0&\cdots&0&a_{n}&b_{n}\end{array}\right|\qquad\left(n\geq 0\right),$

with for $n=0$ the notation

(2.2) $K\binom{\ast}{b_{0}}:=b_{0}.$

Developing with respect to the last line yields for $n\geq 0$

(2.3)
$K\binom{a_{1},\ldots,a_{n+2}}{b_{0},\ldots,b_{n+2}}=b_{n+2}K\binom{a_{1},\ldots,a_{n+1}}{b_{0},\ldots,b_{n+1}}+a_{n+2}K\binom{a_{1},\ldots,a_{n}}{b_{0},\ldots,b_{n}},$

Developing with respect to the first column yields for $n\geq 0$

(2.4)
$K\binom{a_{1},\ldots,a_{n+2}}{b_{0},\ldots,b_{n+2}}=b_{0}K\binom{a_{2},\ldots,a_{n+2}}{b_{1},\ldots,b_{n+2}}+a_{1}K\binom{a_{3},\ldots,a_{n+2}}{b_{2},\ldots,b_{n+2}}.$

Let $A_{n}$ and $B_{n}$ be defined in (1.1), (1.2) and (1.3). It is clear from
(2.1) that

(2.5) $\displaystyle K\binom{\ast}{b_{0}}$ $\displaystyle=b_{0}=A_{0},\qquad
K\binom{a_{1}}{b_{0},b_{1}}=a_{1}+b_{0}b_{1}=A_{1},$ (2.6) $\displaystyle
K\binom{\ast}{b_{1}}$ $\displaystyle=b_{1}=B_{1},\qquad
K\binom{a_{2}}{b_{1},b_{2}}=a_{2}+b_{1}b_{2}=B_{2}.$

Hence an easy induction using (2.3) shows that

(2.7) $\displaystyle A_{n}$
$\displaystyle=K\binom{a_{1},\ldots,a_{n}}{b_{0},\ldots,b_{n}}\qquad\left(n\geq
0\right),$ (2.8) $\displaystyle B_{n}$
$\displaystyle=K\binom{a_{2},\ldots,a_{n}}{b_{1},\ldots,b_{n}}\qquad\left(n\geq
1\right).$

As a first application, we have by following [10, Page 9]:

######  Proposition 1.

Let $n\geq 0.$ Define $A_{n},$ $B_{n},$ $A_{n}^{\prime}$ and $B_{n}^{\prime}$
by

$\displaystyle\frac{A_{n}}{B_{n}}$
$\displaystyle=b_{0}+\frac{a_{1}}{b_{1}}\genfrac{}{}{0.0pt}{}{{}}{+}\frac{a_{2}}{b_{2}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots+}\frac{a_{n}}{b_{n}},$
$\displaystyle\frac{A_{n}^{\prime}}{B_{n}^{\prime}}$
$\displaystyle=b_{n}+\frac{a_{n}}{b_{n-1}}\genfrac{}{}{0.0pt}{}{{}}{+}\frac{a_{n-1}}{b_{n-2}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots+}\frac{a_{1}}{b_{0}}.$

Then for all $n\geq 1,$

(2.9) $A_{n}^{\prime}=A_{n},\quad B_{n}^{\prime}=A_{n-1},\quad
A_{n-1}^{\prime}=B_{n},\quad B_{n-1}^{\prime}=B_{n-1}.$

###### Proof.

We have by (2.7) and (2.8)

$\displaystyle A_{n}^{\prime}$
$\displaystyle=K\binom{a_{n},\ldots,a_{1}}{b_{n},\ldots,b_{0}}=K\binom{a_{1},\ldots,a_{n}}{b_{0},\ldots,b_{n}}=A_{n},$
$\displaystyle B_{n}^{\prime}$
$\displaystyle=K\binom{a_{n-1},\ldots,a_{1}}{b_{n-1},\ldots,b_{0}}=K\binom{a_{1},\ldots,a_{n-1}}{b_{0},\ldots,b_{n-1}}=A_{n-1},$
$\displaystyle A_{n-1}^{\prime}$
$\displaystyle=K\binom{a_{n},\ldots,a_{2}}{b_{n},\ldots,b_{1}}=K\binom{a_{2},\ldots,a_{n}}{b_{1},\ldots,b_{n}}=B_{n},$
$\displaystyle B_{n-1}^{\prime}$
$\displaystyle=K\binom{a_{n-1},\ldots,a_{2}}{b_{n-1},\ldots,b_{1}}=K\binom{a_{2},\ldots,a_{n-1}}{b_{1},\ldots,b_{n-1}}=B_{p-1},$

since a determinant is unchanged by symmetry with respect to its second
diagonal222Up to now, I didn’t know this result… I could find by myself an
elementary proof of it, but I have no reference. Do you know one?. ∎

More generally, define

(2.10)
$\alpha_{k,n}:=b_{k}+\frac{a_{k+1}}{b_{k+1}}\genfrac{}{}{0.0pt}{}{{}}{+\cdots+}\frac{a_{k+n}}{b_{k+n}}=\frac{A_{k,n}}{B_{k,n}}\qquad\left(k\geq
0,n\geq 1\right).$

Then $\alpha_{0,n}=\alpha_{n},$ $A_{0,n}=A_{n},$ $B_{0,n}=B_{n}$ and

(2.11) $A_{k,-1}=1,\quad A_{k,0}=b_{k},\quad B_{k,-1}=0,\quad B_{k,0}=1.$

Moreover, by (2.7) and (2.8)

(2.12) $\displaystyle A_{k,n}$
$\displaystyle=K\binom{a_{k+1},\ldots,a_{k+n}}{b_{k},\ldots,b_{k+n}}\qquad\left(k\geq
0,n\geq 0\right),$ (2.13) $\displaystyle B_{k,n}$
$\displaystyle=K\binom{a_{k+2},\ldots,a_{k+n}}{b_{k+1},\ldots,b_{k+n}}\qquad\left(k\geq
0,n\geq 1\right).$

As $B_{k,0}=A_{k+1,-1}=1,$ this yields immediately

(2.14) $B_{k,n}=A_{k+1,n-1}\qquad\left(k\geq 0,n\geq 0\right).$

Clearly $A_{k,1}=b_{k}A_{k+1,0}+a_{k+1}A_{k+2,-1}$ by (2.5) and (2.11).
Therefore replacing $a_{n}$ and $b_{n}$ by $a_{k+n}$ and $b_{k+n}$ in (2.4)
yields

(2.15) $A_{k,n+2}=b_{k}A_{k+1,n+1}+a_{k+1}A_{k+2,n}\qquad\left(k\geq
0,n\geq-1\right).$

So by using (2.14), (2.6) and (2.11)

(2.16) $B_{k,n+2}=b_{k+1}B_{k+1,n+1}+a_{k+2}B_{k+2,n}\qquad\left(k\geq
0,n\geq-1\right).$

######  Proposition 2.

If $k\geq 0$ and $n\geq 1,$ then

(2.17) $\displaystyle A_{n+k}$
$\displaystyle=A_{n,k}A_{n-1}+a_{n}B_{n,k}A_{n-2},$ (2.18) $\displaystyle
B_{n+k}$ $\displaystyle=A_{n,k}B_{n-1}+a_{n}B_{n,k}B_{n-2},$

###### Proof.

We follow [9, Prop.1]. Let $m=n+k.$ We have to prove that

(2.19) $\displaystyle A_{m}$
$\displaystyle=A_{n,m-n}A_{n-1}+a_{n}B_{n,m-n}A_{n-2}\qquad\left(1\leq n\leq
m\right),$ (2.20) $\displaystyle B_{m}$
$\displaystyle=A_{n,m-n}B_{n-1}+a_{n}B_{n,m-n}B_{n-2}\qquad\left(1\leq n\leq
m\right).$

Let $f(n):=A_{n,m-n}A_{n-1}+a_{n}B_{n,m-n}A_{n-2}.$ Then for $1\leq n<m$

$\displaystyle f(n+1)$
$\displaystyle=A_{n+1,m-n-1}A_{n}+a_{n+1}B_{n+1,m-n-1}A_{n-1}$
$\displaystyle=A_{n+1,m-n-1}(b_{n}A_{n-1}+a_{n}A_{n-2})+a_{n+1}B_{n+1,m-n-1}A_{n-1}$
$\displaystyle=\left(b_{n}A_{n+1,m-n-1}+a_{n+1}B_{n+1,m-n-1}\right)A_{n-1}+a_{n}A_{n+1,m-n-1}A_{n-2}$
$\displaystyle=\left(b_{n}A_{n+1,m-n-1}+a_{n+1}A_{n+2,m-n-2}\right)A_{n-1}+a_{n}B_{n,m-n}A_{n-2}$
$\displaystyle=A_{n,m-n}A_{n-1}+a_{n}B_{n,m-n}A_{n-2}=f(n).$

Hence
$f(n)=f(m)=A_{n,0}A_{n-1}+a_{n}B_{n,0}A_{n-2}=b_{n}A_{n-1}+a_{n}A_{n-2}=A_{n},$
which proves (2.19). The proof of (2.20) is similar: just replace $A_{n},$
$A_{n-1}$ and $A_{n-2}$ by $B_{n},$ $B_{n-1}$ and $B_{n-2}$ respectively. ∎

As $A_{n,0}=b_{n}$ and $B_{n,0}=1$ by (2.11), (1.2) and (1.3) result from
(2.17) and (2.18) respectively by taking $k=0$.

######  Proposition 3.

Let $n\geq 1$ and $k\geq 0.$ Then

(2.21) $A_{n+k}B_{n-1}-A_{n-1}B_{n+k}=\left(-1\right)^{n-1}a_{1}a_{2}\cdots
a_{n}B_{n,k}.$

###### Proof.

For $n\geq 1$ and $k\geq 0$, let $h(n,k):=A_{n+k}B_{n-1}-A_{n-1}B_{n+k}.$Then
by (2.17) and (2.18)

$\displaystyle h(n,k)$
$\displaystyle=\left(A_{n,k}A_{n-1}+a_{n}B_{n,k}A_{n-2}\right)B_{n-1}-A_{n-1}\left(A_{n,k}B_{n-1}+a_{n}B_{n,k}B_{n-2}\right)$
$\displaystyle=a_{n}B_{n,k}\left(A_{n-2}B_{n-1}-A_{n-1}B_{n-2}\right).$

If $n=1,$ then (2.21) is true by (1.1). If $n\geq 2,$ (2.21) results from
(1.5). ∎

From (2.21) we deduce immediately the following generalization of (1.6):

(2.22)
$\frac{A_{n+k}}{B_{n+k}}-\frac{A_{n-1}}{B_{n-1}}=\left(-1\right)^{n-1}a_{1}a_{2}\cdots
a_{n}\frac{B_{n,k}}{B_{n-1}B_{n+k}}\qquad\left(n\geq 1,k\geq 0\right).$

Since $B_{n,0}=1$ by (2.11), (1.5) and (1.6) result from (2.21) and (2.22)
respectively by taking $k=0$. As observed by Perron, for $k=1$ we obtain by
(2.14) and (2.11)

(2.23)
$\frac{A_{n+1}}{B_{n+1}}-\frac{A_{n-1}}{B_{n-1}}=\left(-1\right)^{n-1}a_{1}a_{2}\cdots
a_{n}\frac{b_{n+1}}{B_{n-1}B_{n+1}}\qquad\left(n\geq 1\right).$

This formula can also be deduced directly from (1.5) and (1.3) by writing

$\displaystyle\frac{A_{n+1}}{B_{n+1}}-\frac{A_{n-1}}{B_{n-1}}$
$\displaystyle=\frac{A_{n+1}}{B_{n+1}}-\frac{A_{n}}{B_{n}}+\frac{A_{n}}{B_{n}}-\frac{A_{n-1}}{B_{n-1}}$
$\displaystyle=\left(-1\right)^{n-1}a_{1}a_{2}\cdots
a_{n}\left(\frac{1}{B_{n-1}B_{n}}-\frac{a_{n+1}}{B_{n}B_{n+1}}\right)$
$\displaystyle=\left(-1\right)^{n-1}a_{1}a_{2}\cdots
a_{n}\frac{B_{n+1}-a_{n+1}B_{n-1}}{B_{n-1}B_{n}B_{n+1}}.$

##  3\. Perron’s proof of Tietze theorem

It makes use of two lemmas.

######  Lemma 1.

Let $\alpha$ be the semi-regular continued fraction defined by (1.7). Then

(3.1) $1\leq B_{k,n}\leq B_{k-1,n+1}\leq B_{k+n}\qquad\left(k\geq 1,n\geq
0\right).$

###### Proof.

[10, p.149-150] Replacing $k$ by $k-1$ and $n$ by $n-1$ in (2.16) yields

(3.2)
$B_{k-1,n+1}-B_{k,n}=\left(b_{k}-1\right)B_{k,n}+a_{k+1}B_{k+1,n-1}\qquad\left(k\geq
1,n\geq 0\right).$

We prove first by induction on $n$ that $B_{k,n}\geq 0$ for $k\geq 0$ and
$n\geq-1.$ This is true for $n=-1$ and $n=0$ for all $k\geq 0$ by (2.11).
Assuming that it is true for $n-1$ and $n,$ we get by (1.8)

$B_{k-1,n+1}-B_{k,n}\geq\left(b_{k}-1\right)\left(B_{k,n}-B_{k+1,n-1}\right)\qquad\left(k\geq
1\right),$

which proves that $B_{k,n+1}\geq 0$ for all $k\geq 0.$ Hence we have

(3.3)
$B_{k-1,n+1}-B_{k,n}\geq\left(b_{k}-1\right)\left(B_{k,n}-B_{k+1,n-1}\right)\qquad\left(k\geq
1,n\geq 0\right).$

However $B_{k+n,0}=1$ and $B_{k+n+1,-1}=0$ by (2.11), so that by an easy
induction

$B_{k-1,n+1}\geq B_{k,n}\qquad\left(k\geq 1,n\geq 0\right).$

This proves (3.1) since $B_{k+n,0}=1$ and $B_{0,n+k}=B_{n+k}.$ ∎

######  Lemma 2.

Let $\alpha$ be the semi-regular continued fraction defined by (1.7). Then

$\lim_{n\rightarrow\infty}B_{n}=+\infty.$

###### Proof.

[10, p.150-151] Let $k\geq 1$ be fixed. If $a_{k+1}=1,$ then by (3.2)

(3.4) $B_{k-1,n+1}-B_{k,n}\geq B_{k+1,n-1}\geq 1\qquad\left(n\geq 1\right).$

On the other hand, if $a_{k+1}=-1,$ then by (3.2) and (1.8)

(3.5) $B_{k-1,n+1}-B_{k,n}\geq B_{k,n}-B_{k+1,n-1}\qquad\left(n\geq 0\right).$

So in both cases

(3.6) $B_{k-1,n+1}-B_{k,n}\geq\min\left(1,B_{k,n}-B_{k+1,n-1}\right).$

Consequently, if $a_{k+1}=1,$ by an easy induction we see that

$B_{k-m,n+m}-B_{k-m+1,n+m-1}\geq 1\qquad\left(1\leq m\leq k\right).$

Summing for $m=1$ to $k$ yields $B_{0,n+k}-B_{k,n}\geq k$ and therefore

(3.7) $a_{k+1}=1\quad\Rightarrow\quad B_{k+n}\geq k+1\qquad\left(k\geq 1,n\geq
1\right).$

If $a_{k+1}=-1,$ then $B_{k-1,1}-B_{k,0}\geq 1$ by taking $n=0$ in (3.5).
Therefore by (3.6)

$B_{k-m,m}-B_{k-m+1,m-1}\geq 1\qquad\left(1\leq m\leq k\right).$

Hence summing again for $m=1$ to $k$ yields $B_{0,k}-B_{k,0}\geq k,$ so that

(3.8) $a_{k+1}=-1\quad\Rightarrow\quad B_{k}\geq k+1\qquad\left(k\geq
1\right).$

Now if $a_{k}=-1$ for all large $n,$ then clearly
$\lim_{k\rightarrow\infty}B_{k}=+\infty$ by (3.8). On the other hand, if there
exist infinitely many $k$ such that $a_{k}=1,$ then
$\lim_{k\rightarrow\infty}B_{k}=+\infty$ by (3.7) and the proof of Lemma 2 is
complete.∎

Now we prove Tietze theorem by following [10, p.151]. By (2.22) and (3.1), we
have

$\left|\frac{A_{n+k}}{B_{n+k}}-\frac{A_{n-1}}{B_{n-1}}\right|\leq\frac{B_{n,k}}{B_{n-1}B_{n+k}}\leq\frac{1}{B_{n-1}}\qquad\left(n\geq
1,k\geq 0\right).$

As $\lim_{n\rightarrow\infty}B_{n}=\infty$ by Lemma 2, $A_{n}/B_{n}$ is a
Cauchy sequence, which proves Tietze theorem.

######  Remark 1.

Two different alternative proofs of Lemma 2, not using continuants, can be
found in [3, Lem.3] and in [9, Th.1(i)].

##  4\. Proof of Theorem 2

We will need the following lemma, which is closely connected to the method of
power iteration in numerical analysis.

######  Lemma 3.

Let $\left(u_{0},v_{0}\right)\in\mathbb{C}^{2}\smallsetminus\left(0,0\right)$
and

$M=\left(\begin{array}[c]{cc}a&c\\\
b&d\end{array}\right),\quad\left(a,b,c,d\right)\in\mathbb{C}^{4},\quad b\neq
0,\quad ad-bc\neq 0.$

Let $\lambda_{1}$ and $\lambda_{2}$ with
$\left|\lambda_{1}\right|\geq\left|\lambda_{2}\right|$ be the eigenvalues of
$M,$ and let $(u_{n},v_{n})$ be defined recursively by

$\left(\begin{array}[c]{c}u_{n+1}\\\
v_{n+1}\end{array}\right)=M\left(\begin{array}[c]{c}u_{n}\\\
v_{n}\end{array}\right)\quad\left(n\geq 0\right).$

Define $\left(x_{1},x_{2}\right)$ by $\lambda_{1}=bx_{1}+d$ and
$\lambda_{2}=bx_{2}+d$. Then:

(i) If $\left|\lambda_{1}\right|>\left|\lambda_{2}\right|,$ $v_{n}\neq 0$ for
all large $n$ and

(4.1) $\displaystyle\lim_{n\rightarrow\infty}\frac{u_{n}}{v_{n}}$
$\displaystyle=x_{1}\quad\text{if\quad}u_{0}-v_{0}x_{2}\neq 0,$ (4.2)
$\displaystyle\lim_{n\rightarrow\infty}\frac{u_{n}}{v_{n}}$
$\displaystyle=x_{2}\quad\text{if\quad}u_{0}-v_{0}x_{2}=0.$

(ii) If $\left|\lambda_{1}\right|=\left|\lambda_{2}\right|$ and
$\lambda_{1}\neq\lambda_{2},$ $u_{n}/v_{n}$ is divergent, except if

$u_{0}-v_{0}x_{1}=0\quad or\quad u_{0}-v_{0}x_{2}=0.$

(iii) If $\lambda_{1}=\lambda_{2},$ $v_{n}\neq 0$ for all large $n$ and
$\lim_{n\rightarrow\infty}u_{n}/v_{n}=x_{1}.$

###### Proof.

We note first that $\lambda_{1}\neq 0$ and $\lambda_{2}\neq 0$ since $\det
M=ad-bc\neq 0.$ Let $f$ be the linear transformation of $\mathbb{C}^{2}$ whose
matrix in the standard basis is $M.$ Let $e_{1},$ $e_{2},$ and $w_{n}$
$\left(n\geq 0\right)$ be defined by

$e_{1}=\left(\begin{array}[c]{c}x_{1}\\\ 1\end{array}\right),\quad
e_{2}=\left(\begin{array}[c]{c}x_{2}\\\ 1\end{array}\right),\quad
w_{n}=\left(\begin{array}[c]{c}u_{n}\\\ v_{n}\end{array}\right).$

It is easy to check that $Me_{1}=\lambda_{1}e_{1}$ and
$Me_{2}=\lambda_{2}e_{2}.$ Hence $e_{1}$ and $e_{2}$ are eigenvectors of $M$
and $f$ and they form a basis of $\mathbb{C}^{2}$ if
$\lambda_{1}\neq\lambda_{2}$. In this case, define $\mu_{1}$ and $\mu_{2}$ by

(4.3) $w_{0}=\mu_{1}e_{1}+\mu_{2}e_{2}.$

We observe that

$\mu_{1}=0\Leftrightarrow\left|\begin{array}[c]{cc}u_{0}&x_{2}\\\
v_{0}&1\end{array}\right|=0\Leftrightarrow u_{0}-v_{0}x_{2}=0$

and similarly $\mu_{2}=0\Leftrightarrow u_{0}-v_{0}x_{1}=0.$ It results from
(4.3) that

(4.4)
$w_{n}=f^{n}\left(w_{0}\right)=\mu_{1}\lambda_{1}^{n}e_{1}+\mu_{2}\lambda_{2}^{n}e_{2}\quad\left(n\geq
0\right).$

Assume that $\left|\lambda_{1}\right|>\left|\lambda_{2}\right|.$ If
$\mu_{1}\neq 0,$ we can write

$w_{n}=\mu_{1}\lambda_{1}^{n}\left(e_{1}+\frac{\mu_{2}}{\mu_{1}}\left(\frac{\lambda_{2}}{\lambda_{1}}\right)^{n}e_{2}\right),$

which proves that $v_{n}\neq 0$ for all large $n$ and (4.1) holds. If
$\mu_{1}=0,$ then (4.2) holds by (4.4), which proves (i). Now assume that
$\left|\lambda_{1}\right|=\left|\lambda_{2}\right|$ and
$\lambda_{1}\neq\lambda_{2}.$ Then $\lambda_{2}=\lambda_{1}e^{i\theta}$ for
some $\theta\in\left]0,2\pi\right[$ and by (4.4) we see that

(4.5)
$w_{n}=\lambda_{1}^{n}\left(\mu_{1}e_{1}+\mu_{2}e^{in\theta}e_{2}\right)\quad\left(n\geq
0\right).$

So if $\mu_{1}=0$ or $\mu_{2}=0$ the sequence $u_{n}/v_{n}$ is well-defined
and convergent. If $\mu_{1}\neq 0$ and $\mu_{2}\neq 0,$ then either $v_{n}=0$
for infinitely many $n$ and $u_{n}/v_{n}$ is not defined, or

$\frac{u_{n}}{v_{n}}=\frac{\mu_{1}x_{1}+\mu_{2}e^{in\theta}x_{2}}{\mu_{1}+\mu_{2}e^{in\theta}}=x_{1}+\left(x_{2}-x_{1}\right)\frac{\mu_{2}}{\mu_{1}e^{-in\theta}+\mu_{2}},$

which is divergent since $x_{2}\neq x_{1},$ $\mu_{1}\neq 0,$ $\mu_{2}\neq 0$
and $\theta\in\left]0,2\pi\right[.$ This proves (ii). Finally, if
$\lambda_{1}=\lambda_{2},$ there exists a basis
$\left(e_{1},e_{2}^{\prime}\right)$ of $\mathbb{C}^{2}$ such that
$f(e_{2}^{\prime})=e_{1}+\lambda_{1}e_{2}^{\prime}$ and so

(4.6) $w_{0}=\mu_{1}e_{1}+\mu_{2}e_{2}^{\prime}$

for some $(\mu_{1},\mu_{2})\in\mathbb{C}^{2}$. As
$f^{n}(e_{2}^{\prime})=n\lambda_{1}^{n-1}e_{1}+\lambda_{1}^{n}e_{2}^{\prime},$
we get

$w_{n}=f^{n}(w_{0})=\lambda_{1}^{n-1}\left(\left(\mu_{1}\lambda_{1}+n\mu_{2}\right)e_{1}+\lambda_{1}\mu_{2}e_{2}^{\prime}\right).$

If $\mu_{2}=0,$ we see that $u_{n}/v_{n}=x_{1}.$ If $\mu_{2}\neq 0,$

$\lim_{n\rightarrow\infty}\frac{u_{n}}{v_{n}}=\lim_{n\rightarrow\infty}\frac{\left(\mu_{1}\lambda_{1}+n\mu_{2}\right)x_{1}+\lambda_{1}\mu_{2}x_{2}}{\left(\mu_{1}\lambda_{1}+n\mu_{2}\right)+\lambda_{1}\mu_{2}}=x_{1}$

and (iii) is proved.∎

Now we prove Theorem 2. Let $q\in\left\\{0,\ldots,p-1\right\\}.$ Replacing $n$
by $p$ and $k$ by $np+q$ in (2.17) and (2.18) yields

$\displaystyle A_{\left(n+1\right)p+q}$
$\displaystyle=A_{p-1}A_{p,np+q}+a_{p}A_{p-2}B_{p,np+q},$ $\displaystyle
B_{\left(n+1\right)p+q}$
$\displaystyle=B_{p-1}A_{p,np+q}+a_{p}B_{p-2}B_{p,np+q}.$

However by (2.12) and (2.7) we have since $a_{n}$ and $b_{n}$ are p-periodic

$A_{p,np+q}=K\binom{a_{p+1},\ldots,a_{\left(n+1\right)p+q}}{b_{p},\ldots,b_{\left(n+1\right)p+q}}=K\binom{a_{1},\ldots,a_{np+q}}{b_{0},\ldots,b_{np+q}}=A_{np+q}.$

Similarly, by (2.13) and (2.8), we see that $B_{p,np+q}=B_{np+q}.$ Hence for
$0\leq q\leq p-1$

(4.7) $\displaystyle A_{\left(n+1\right)p+q}$
$\displaystyle=A_{p-1}A_{np+q}+a_{p}A_{p-2}B_{np+q},$ (4.8) $\displaystyle
B_{\left(n+1\right)p+q}$ $\displaystyle=B_{p-1}A_{np+q}+a_{p}B_{p-2}B_{np+q}.$

Assume that $B_{p-1}=0.$ Then by (4.8) with $q=p-1$ we have

$B_{\left(n+1\right)p+p-1}=a_{p}B_{p-2}B_{np+p-1}=\left(a_{p}B_{p-2}\right)^{n+1}B_{p-1}=0,$

so that the sequence $B_{n}$ vanishes infinitely often and therefore $\alpha$
is divergent.

Assume that $B_{p-1}\neq 0.$ Then (4.7) and (4.8) can be written as

$\left(\begin{array}[c]{c}A_{\left(n+1\right)p+q}\\\
B_{\left(n+1\right)p+q}\end{array}\right)=\left(\begin{array}[c]{cc}A_{p-1}&a_{p}A_{p-2}\\\
B_{p-1}&a_{p}B_{p-2}\end{array}\right)\left(\begin{array}[c]{c}A_{np+q}\\\
B_{np+q}\end{array}\right)=M\left(\begin{array}[c]{c}A_{np+q}\\\
B_{np+q}\end{array}\right).$

We observe that $\det M=\left(-1\right)^{p}a_{1}\cdots a_{p}\neq 0,$ so that
we can apply Lemma 3 with $u_{n}=u_{n,q}=A_{np+q}$ and
$v_{n}=v_{n,q}=B_{np+q}$ for every $q\in\left\\{0,\ldots,p-1\right\\},$ in
such a way that

$u_{0,q}-v_{0,q}x_{1}=A_{q}-x_{1}B_{q},\quad
u_{0,q}-v_{0,q}x_{2}=A_{q}-x_{2}B_{q}.$

As $\lambda_{1}+\lambda_{2}=A_{p-1}+a_{p}B_{p-2},$ for $q=p-1$ we see that

(4.9) $\displaystyle u_{0,p-1}-v_{0,p-1}x_{1}$
$\displaystyle=A_{p-1}-x_{1}B_{p-1}=A_{p-1}-\left(\lambda_{1}-a_{p}B_{p-2}\right)=\lambda_{2}\neq
0,$ (4.10) $\displaystyle u_{0,p-1}-v_{0,p-1}x_{2}$
$\displaystyle=A_{p-1}-x_{2}B_{p-1}=A_{p-1}-\left(\lambda_{2}-a_{p}B_{p-2}\right)=\lambda_{1}\neq
0.$

We distinguish three cases.

Case 1. $\lambda_{1}=\lambda_{2}.$ By Lemma 3 (iii) we see that
$\lim_{n\rightarrow\infty}A_{np+q}/B_{np+q}=x_{1}$ for every
$q\in\left\\{0,\ldots,p-1\right\\},$ so that $\alpha$ is convergent and
$\alpha=x_{1}.\vskip 3.0pt plus 1.0pt minus 1.0pt$

Case 2. $\left|\lambda_{1}\right|=\left|\lambda_{2}\right|$ and
$\lambda_{1}\neq\lambda_{2}.$ Then the sequence $A_{np+p-1}/B_{np+p-1}$ is
divergent by Lemma 3 (ii), (4.9) and (4.10). Therefore $\alpha$ is divergent.

Case 3. $\left|\lambda_{1}\right|>\left|\lambda_{2}\right|.$ Then the sequence
$A_{np+p-1}/B_{np+p-1}$ converges to $x_{1}$ by Lemma 3 (i) and (4.10).
Moreover for $0\leq q\leq p-2$ the sequence $A_{np+p-1}/B_{np+p-1}$ converges
to $x_{1}$ if $A_{q}-x_{2}B_{q}\neq 0$ and to $x_{2}$ if $A_{q}-x_{2}B_{q}=0.$
As $x_{1}\neq x_{2}$ since $\lambda_{1}\neq\lambda_{2},$ $\alpha$ is
convergent if and only if $A_{q}-x_{2}B_{q}\neq 0$ for all $0\leq q\leq p-2,$
and the proof of Theorem 2 is complete.

######  Remark 2.

In Case 3, the sequence $A_{n}/B_{n}$ has two different limit points
$x_{1}\neq x_{2}$ when $A_{q}-x_{2}B_{q}=0$ for some
$q\in\left\\{0,\ldots,p-2\right\\}$ ([14],[11, Satz 2.39],[5, Theorem 3.3
(B)]). This is known as Thiele’s oscillations.

##  5\. Proofs of Theorem 3 and Corollary 1

First we prove Theorem 3. By Proposition 1, we have

$A_{p}^{\prime}=A_{p},\quad B_{p}^{\prime}=A_{p-1},\quad
A_{p-1}^{\prime}=B_{p},\quad B_{p-1}^{\prime}=B_{p-1}.$

Hence the matrix $M^{\prime}$ associated to $\alpha^{\prime}$ in Theorem 2 is

$M^{\prime}:=\left(\begin{array}[c]{cc}A_{p-1}^{\prime}&A_{p}^{\prime}-b_{0}^{\prime}A_{p-1}^{\prime}\\\
B_{p-1}^{\prime}&B_{p}^{\prime}-b_{0}^{\prime}B_{p-1}^{\prime}\end{array}\right)=\left(\begin{array}[c]{cc}B_{p}&A_{p}-b_{0}B_{p}\\\
B_{p-1}&A_{p-1}-b_{0}B_{p-1}\end{array}\right).$

Its trace is
$\operatorname*{tr}M^{\prime}=B_{p}+A_{p-1}-b_{0}B_{p-1}=\operatorname*{tr}M$
and its determinant is

$\det M^{\prime}=A_{p-1}B_{p}-A_{p}B_{p-1}=\det M.$

So the eigenvalues of $M^{\prime}$ are exactly
$\lambda_{1}^{\prime}=\lambda_{1}$ and $\lambda_{2}^{\prime}=\lambda_{2}.$ Now
let $x_{1}^{\prime}$ and $x_{2}^{\prime}$ be defined by

$\lambda_{1}^{\prime}=\lambda_{1}=x_{1}^{\prime}B_{p-1}^{\prime}+a_{p}^{\prime}B_{p-2}^{\prime},\quad\lambda_{2}^{\prime}=\lambda_{2}=x_{2}^{\prime}B_{p-1}^{\prime}+a_{p}^{\prime}B_{p-2}^{\prime}.$

Since $b_{p}^{\prime}=b_{p}=b_{0,}$ we see by (4.9) and (4.10) that

$\displaystyle\left(b_{0}-x_{1}\right)B_{p-1}^{\prime}+a_{p}^{\prime}B_{p-2}^{\prime}$
$\displaystyle=\left(b_{0}-x_{1}\right)B_{p-1}^{\prime}+B_{p}^{\prime}-b_{0}B_{p-1}^{\prime}=\lambda_{2},$
$\displaystyle\left(b_{0}-x_{2}\right)B_{p-1}^{\prime}+a_{p}^{\prime}B_{p-2}^{\prime}$
$\displaystyle=\left(b_{0}-x_{2}\right)B_{p-1}^{\prime}+B_{p}^{\prime}-b_{0}B_{p-1}^{\prime}=\lambda_{1}.$

Therefore $x_{1}^{\prime}=b_{0}-x_{2}$ and $x_{2}^{\prime}=b_{0}-x_{1},$ so
that Theorem 3 follows immediately from Theorem 2.

Now we prove Corollary 1. As $\alpha=x_{1}$ is irrational, (1.13) cannot hold
and the reverse continued fraction $\alpha^{\prime}$ of $\alpha$ converges to
$b_{0}-x_{2}$ by Theorem 3. Moreover $\lambda_{1}$ is irrational by (1.10).
Hence $\lambda_{1}$ is quadratic irrational and $\lambda_{2}$ is the conjugate
of $\lambda_{1}$ since the characteristic polynomial of $M$ has integer
coefficients. Consequently $x_{2}$ is the conjugate of $x_{1}=\alpha$ by
(1.9), which yields $\alpha^{\prime}=b_{0}-x_{2}=b_{0}-\alpha^{\ast}.$
Corollary 1 is proved.

## References

  * [1] D. Duverney, Number theory, World Scientific, 2010.
  * [2] D. Duverney and I. Shiokawa, Irrationality exponents of semi-regular continued fractions, Arxiv:2202.10738v1 (22 Feb 2022).
  * [3] D. Duverney and I. Shiokawa, A new proof of a theorem of Tietze, Arxiv:2203.05296v1 (10 Mar 2022).
  * [4] É. Galois, Démonstration d’un théorème sur les fractions continues périodiques, Ann. Math. Pures. Appl. 19 (1828-1829), 294-301.
  * [5] W. B. Jones and W. J. Thron, Continued fractions : analytic theory and applications, Encyclopedia of Math. and its applications 11, Addison-Wesley, 1980.
  * [6] A. F. Möbius, Beiträge zu der Lehre von den Kettenbrüchen nebst einem Anhange dioptrischen Inhalts, J. Reine Angew. Math. 6 (1830), 215-243.
  * [7] T. Muir, A new special class of determinants, Proc. Ed. 8 (1874), 229-236.
  * [8] V. Nachreiner, Beziehungen zwischen Determinanten und Kettenbrüchen, Preisschrift München (1872).
  * [9] A. Offer, Continuants and semi-regular continued fractions (2008), http://www.numbertheory.org/PDFS/continuant.pdf.
  * [10] O. Perron, Die Lehre von den Kettenbrüchen, Dritte Auflage, Band I, Teubner, 1954 (reprint Springer 1977).
  * [11] O. Perron, Die Lehre von den Kettenbrüchen, Dritte Auflage, Band II, Teubner, 1957 (reprint Springer 1977).
  * [12] W. Spottiswoode, Elementary theorems relating to determinants, Crelle J. 51 (1856), 209-271.
  * [13] O. Stolz, Vorlesungen über allgemeine Arithmetik, Leipzig, 1886.
  * [14] T. L. Thiele, Bemaerkninger om periodiske Kjaedebrokers Konvergens, Tidsskrift for Mathematik (4) 3 (1879).
  * [15] H. Tietze, Über Kriterien für Konvergenz und Irrationalität unendlichen Kettenbrüche, Math. Ann. 70 (1911), 236-265.
  * [16] E. Zurl, Theorie der reduziert-regelmässigen Kettenbrüche, Math. Ann. 110, 679-717 (1935).


# Algorithms for Approximate Minimization of the Difference Between Submodular
Functions, with Applications††thanks: A shorter version of this appeared in
Proc. Uncertainty in Artificial Intelligence (UAI), Catalina Islands, 2012

Rishabh Iyer
Dept. of Electrical Engineering
University of Washington
Seattle, WA-98175, USA  Jeff Bilmes
Dept. of Electrical Engineering
University of Washington
Seattle, WA-98175, USA

###### Abstract

We extend the work of Narasimhan and Bilmes [32] for minimizing set functions
representable as a difference between submodular functions. Similar to [32],
our new algorithms are guaranteed to monotonically reduce the objective
function at every step. We empirically and theoretically show that the per-
iteration cost of our algorithms is much less than [32], and our algorithms
can be used to efficiently minimize a difference between submodular functions
under various combinatorial constraints, a problem not previously addressed.
We provide computational bounds and a hardness result on the multiplicative
inapproximability of minimizing the difference between submodular functions.
We show, however, that it is possible to give worst-case additive bounds by
providing a polynomial time computable lower-bound on the minima. Finally we
show how a number of machine learning problems can be modeled as minimizing
the difference between submodular functions. We experimentally show the
validity of our algorithms by testing them on the problem of feature selection
with submodular cost features.

##  1 Introduction

Discrete optimization is important to many areas of machine learning and
recently an ever growing number of problems have been shown to be expressible
as submodular function minimization or maximization (e.g., [21, 25, 27, 30,
29, 31]). The class of submodular functions is indeed special since submodular
function minimization is known to be polynomial time, while submodular
maximization, although NP complete, admits constant factor approximation
algorithms. Let $V=\\{1,2,\cdots,n\\}$ refer a ground set, then
$f:2^{V}\rightarrow\mathbb{R}$ is said to be submodular if for sets
$S,T\subseteq V$, $f(S)+f(T)\geq f(S\cup T)+f(S\cap T)$ (see [11] for details
on submodular, supermodular, and modular functions). Submodular functions have
a diminishing returns property, wherein the gain of an element in the context
of bigger set is lesser than the gain of that element in the context of a
smaller subset. This property occurs naturally in many applications in machine
learning, computer vision, economics, operations research, etc.

In this paper, we address the following problem. Given two submodular
functions $f$ and $g$, and define $v(X)\triangleq f(X)-g(X)$, solve the
following optimization problem:

$\min_{X\subseteq V}[f(X)-g(X)]\equiv\min_{X\subseteq V}[v(X)].$ (1)

A number of machine learning problems involve minimization over a difference
between submodular functions. The following are some examples:

  * •

Sensor placement with submodular costs:  The problem of choosing sensor
locations $A$ from a given set of possible locations $V$ can be modeled [25,
26] by maximizing the mutual information between the chosen variables $A$ and
the unchosen set $V\backslash A$ (i.e., $f(A)=I(X_{A};X_{V\backslash A})$).
Alternatively, we may wish to maximize the mutual information between a set of
chosen sensors $X_{A}$ and a fixed quantity of interest $C$ (i.e.,
$f(A)=I(X_{A};C)$) under the assumption that the set of features $X_{A}$ are
conditionally independent given $C$ [25]. These objectives are submodular and
thus the problem becomes maximizing a submodular function subject to a
cardinality constraint. Often, however, there are costs $c(A)$ associated with
the locations that naturally have a diminishing returns property. For example,
there is typically a discount when purchasing sensors in bulk. Moreover, there
may be diminished cost for placing a sensor in a particular location given
placement in certain other locations (e.g., the additional equipment needed to
install a sensor in, say, a precarious environment could be re-used for
multiple sensor installations in like environments). Hence, along with
maximizing mutual information, we also want to simultaneously minimize the
cost and this problem can be addressed by minimizing the difference between
submodular functions $f(A)-\lambda c(A)$ for tradeoff parameter $\lambda$.

  * •

Discriminatively structured graphical models and neural computation:  An
application suggested in [32] and the initial motivation for this problem is
to optimize the EAR criterion to produce a discriminatively structured
graphical model. EAR is basically a difference between two mutual information
functions (i.e., a difference between submodular functions). [32] shows how
classifiers based on discriminative structure using EAR can significantly
outperform classifiers based on generative graphical models. Note also that
the EAR measure is the same as “synergy” in a neural code [3], widely used in
neuroscience.

  * •

Feature selection:  Given a set of features $X_{1},X_{2},\cdots,X_{|V|}$, the
feature selection problem is to find a small subset of features $X_{A}$ that
work well when used in a pattern classifier. This problem can be modeled as
maximizing the mutual information $I(X_{A};C)$ where $C$ is the class. Note
that $I(X_{A};C)=H(X_{A})-H(X_{A}|C)$ is always a difference between
submodular functions. Under the naïve Bayes model, this function is submodular
[25]. It is not submodular under general classifier models such as support
vector machines (SVMs) or neural networks. Certain features, moreover, might
be cheaper to use given that others are already being computed. For example,
if a subset $S_{i}\subseteq V$ of the features for a particular information
source $i$ are spectral in nature, then once a particular $v\in S_{i}$ is
chosen, the remaining features $S_{i}\setminus\\{v\\}$ may be relatively
inexpensive to compute, due to grouped computational strategies such as the
fast Fourier transform. Therefore, it might be more appropriate to use a
submodular cost model $c(A)$. One such cost model might be
$c(A)=\sum_{i}\sqrt{m(A\cap S_{i})}$ where $m(j)$ would be the cost of
computing feature $j$. Another might be $c(A)=\sum_{i}c_{i}\min(|A\cap
S_{i}|,1)$ where $c_{i}$ is the cost of source $i$. Both offer diminishing
cost for choosing features from the same information source. Such a cost model
could be useful even under the naïve Bayes model, where $I(X_{A};C)$ is
submodular. Feature selection becomes a problem of maximizing
$I(X_{A};C)-\lambda c(A)=H(X_{A})-[H(X_{A}|C)+\lambda c(A)]$, the difference
between two submodular functions.

  * •

Probabilistic Inference: A typical instance of probabilistic inference is the
following: We are given a distribution $p(x)\propto\exp(-v(x))$ where
$x\in\\{0,1\\}^{n}$ and $v$ is a pseudo-Boolean function [2]. It is desirable
to compute $\operatorname*{argmax}_{x\in\\{0,1\\}^{n}}p(x)$ which means
minimizing $v(x)$ over $x$, the most-probable explanation (MPE) problem [35].
If $p$ factors with respect to a graphical model of tree-width $k$, then
$v(x)=\sum_{i}v_{i}(x_{C})$ where $C_{i}$ is a bundle of indices such that
$|C|\leq k+1$ and the sets $\mathcal{C}=\\{C_{i}\\}_{i}$ form a junction tree,
and it might be possible to solve inference using dynamic programming. If $k$
is large and/or if hypertree factorization does not hold, then approximate
inference is typically used [40]. On the other hand, defining
$x(X)=\\{x\in\\{0,1\\}^{n}:x_{i}=1\text{ whenever }i\in X\\}$, if the set
function $\bar{v}(X)=v(x(X))$ is submodular, then even if $p$ has large tree-
width, the MPE problem can be solved exactly in polynomial time [18]. This, in
fact, is the basis behind inference in many computer vision models where $v$
is often not only submodular but also has limited sized $|C_{i}|$. For
example, for submodular $v$ and if $|C_{i}|\leq 2$ then graph-cuts can solve
the MPE problem extremely rapidly [24] and even some cases with $v$ non-
submodular [23]. An important challenge is to consider non-submodular $v$ that
can be minimized efficiently and for which there are approximation guarantees,
a problem recently addressed in [19]. On the other hand, if $v$ can be
expressed as a difference between two submodular functions (which it can, see
Lemma 3.1), or if such a decomposition can be computed (which it sometimes
can, see Lemma 3.2), then a procedure to minimize the difference between two
submodular functions offers new ways to solve probabilistic inference. As an
example, a large class of rich higher potentials can be expressed as [13]:

$\displaystyle f(x)=\sum_{C\in\mathcal{C}}w_{C}\prod_{i\in C}x_{i}$ (2)

$\mathcal{C}$ here stands for a set of sets, possibly with higher-order terms
(i.e there exist $C\in\mathcal{C}:|C|>2$). If $w_{C}\leq 0,\forall
C\in\mathcal{C}$, then $f$ is submodular. If $|\mathcal{C}|$ is not large (say
polynomial in $n$), we can efficiently find a decomposition into submodular
components (which will contain the sets $C\in\mathcal{C}:w_{C}\leq 0$) and the
supermodular terms (which contain sets $C\in\mathcal{C}:w_{C}\geq 0$). These
can potentially represent a rich class of potential functions for a number of
applications, particularly in vision.

We note that given a solution to Equation 1, we can also minimize the
difference between two supermodular functions $\min((-g)-(-f))$, maximize the
difference between two submodular functions $\max(-v)=\max(g-f)$, and maximize
the difference between two supermodular functions $\max(-v)=\max((-f)-(-g))$.

Previously, Narasimhan and Bilmes [32] proposed an algorithm inspired by the
convex-concave procedure [41] to address Equation (1). This algorithm
iteratively minimizes a submodular function by replacing the second submodular
function $g$ by it’s modular lower bound. They also show that any set function
can be expressed as a difference between two submodular functions and hence
every set function optimization problem can be reduced to minimizing a
difference between submodular functions. They show that this process converges
to a local minima, however the convergence rate is left as an open question.

In this paper, we first describe tight modular bounds on submodular functions
in Section 2, including lower bounds based on points in the base polytope as
used in [32], and recent upper bounds first described in a result in [17]. In
section 3 Submodular-Supermodular Procedure, we describe the submodular-
supermodular procedure proposed in [32]. We further provide a constructive
procedure for finding the submodular functions $f$ and $g$ for any arbitrary
set function $v$. Although our construction is NP hard in general, we show how
for certain classes of set functions $v$, it is possible to find the
decompositions $f$ and $g$ in polynomial time. In Section 4, we propose two
new algorithms both of which are guaranteed to monotonically reduce the
objective at every iteration and which converge to a local minima. Further we
note that the per-iteration cost of our algorithms is in general much less
than [32], and empirically verify that our algorithms are orders of magnitude
faster on real data. We show that, unlike in [32], our algorithms can be
extended to easily optimize equation (1) under cardinality, knapsack, and
matroid constraints. Moreover, one of our algorithms can actually handle
complex combinatorial constraints, such as spanning trees, matchings, cuts,
etc. Further in Section 5, we give a hardness result that there does not exist
any polynomial time algorithm with any polynomial time multiplicative
approximation guarantees unless P=NP, even when it is easy to find or when we
are given the decomposition $f$ and $g$, thus justifying the need for
heuristic methods to solve this problem. We show, however, that it is possible
to get additive bounds by showing polynomial time computable upper and lower
bound on the optima. We also provide computational bounds for all our
algorithms (including the submodular-supermodular procedure), a problem left
open in [32].

Finally we perform a number of experiments on the feature selection problem
under various cost models, and show how our algorithms used to maximize the
mutual information perform better than greedy selection (which would be near
optimal under the naïve Bayes assumptions) and with less cost.

##  2 Modular Upper and Lower bounds

The Taylor series approximation of a convex function provides a natural way of
providing lower bounds on such a function. In particular the first order
Taylor series approximation of a convex function is a lower bound on the
function, and is linear in $x$ for a given $y$ and hence given a convex
function $\phi$, we have:

$\phi(x)\geq\phi(y)+\langle\nabla\phi(y),x-y\rangle.$ (3)

Surprisingly, any submodular function has both a tight lower [8] and upper
bound [17], unlike strict convexity where there is only a tight first order
lower bound.

###  2.1 Modular Lower Bounds

Recall that for submodular function $f$, the submodular polymatroid, base
polytope and the sub-differential with respect to a set $Y$ [11] are
respectively:

$\displaystyle\mathcal{P}_{f}=\\{x:x(S)\leq f(S),\forall S\subseteq V\\}$ (4)
$\displaystyle\mathcal{B}_{f}=\mathcal{P}_{f}\cap\\{x:x(V)=f(V)\\}$ (5)
$\displaystyle\\!\\!\\!\\!\partial f(Y)=\\{y\in\mathbb{R}^{V}:\forall
X\subseteq V,f(Y)-y(Y)\leq f(X)-y(X)\\}$

The extreme points of this sub-differential are easy to find and characterize,
and can be obtained from a greedy algorithm ([8, 11]) as follows:

######  Theorem 2.1.

([11], Theorem 6.11) A point $y$ is an extreme point of $\partial f(Y)$, iff
there exists a chain $\emptyset=S_{0}\subset S_{1}\subset\cdots\subset S_{n}$
with $Y=S_{j}$ for some $j$, such that $y(S_{i}\setminus
S_{i-1})=y(S_{i})-y(S_{i-1})=f(S_{i})-f(S_{i-1})$.

Let $\sigma$ be a permutation of $V$ and define
$S_{i}^{\sigma}=\\{\sigma(1),\sigma(2),\dots,\sigma(i)\\}$ as $\sigma$’s chain
containing $Y$, meaning $S_{|Y|}^{\sigma}=Y$ (we say that $\sigma$’s chain
contains $Y$). Then we can define a sub-gradient $h^{f}_{Y}$ corresponding to
$f$ as:

$h^{f}_{Y,\sigma}(\sigma(i))=\begin{cases}f(S_{1}^{\sigma})&\text{ if }i=1\\\
f(S_{i}^{\sigma})-f(S_{i-1}^{\sigma})&\text{ otherwise }\end{cases}.$

We get a modular lower bound of $f$ as follows:

$\displaystyle h^{f}_{Y,\sigma}(X)\leq f(X),\forall X\subseteq V,\text{ and
}\forall i,h^{f}_{Y,\sigma}(S_{i}^{\sigma})=f(S_{i}^{\sigma}),$

which is parameterized by a set $Y$ and a permutation $\sigma$. Note
$h(X)=\sum_{i\in X}h(i)$, and $h^{f}_{Y,\sigma}(Y)=f(Y)$. Observe the
similarity to convex functions, where a linear lower bound is parameterized by
a vector $y$.

###  2.2 Modular Upper Bounds

For $f$ submodular, [33] established the following:

$\displaystyle f(Y)\leq f(X)-\sum_{j\in X\backslash Y}f(j|X\backslash
j)+\sum_{j\in Y\backslash X}f(j|X\cap Y),$ $\displaystyle f(Y)\leq
f(X)-\sum_{j\in X\backslash Y}f(j|(X\\!\cup\\!Y)\backslash j)+\sum_{j\in
Y\backslash X}f(j|X)$

Note that $f(A|B)\triangleq f(A\cup B)-f(B)$ is the gain of adding $A$ in the
context of $B$. These upper bounds in fact characterize submodular functions,
in that a function $f$ is a submodular function iff it follows either of the
above bounds. Using the above, two tight modular upper bounds ([17]) can be
defined as follows:

$\displaystyle\\!\\!\\!f(Y)\leq m^{f}_{X,1}(Y)\triangleq
f(X)-\\!\\!\\!\sum_{j\in X\backslash Y}f(j|X\backslash j)+\sum_{j\in
Y\backslash X}f(j|\emptyset),$ $\displaystyle\\!\\!\\!f(Y)\leq
m^{f}_{X,2}(Y)\triangleq f(X)-\\!\\!\\!\sum_{j\in X\backslash
Y}f(j|V\backslash j)+\sum_{j\in Y\backslash X}f(j|X).$

Hence, this yields two tight (at set $X$) modular upper bounds
$m^{f}_{X,1},m^{f}_{X,2}$ for any submodular function $f$. For briefness, when
referring either one we use $m^{f}_{X}$.

## 3 Submodular-Supermodular Procedure

We now review the submodular-supermodular procedure [32] to minimize functions
expressible as a difference between submodular functions (henceforth called DS
functions). Interestingly, any set function can be expressed as a DS function
using suitable submodular functions as shown below. The result was first shown
in [32] using the Lovász extension. We here give a new combinatorial proof,
which avoids Hessians of polyhedral convex functions and which provides a way
of constructing (a non-unique) pair of submodular functions $f$ and $g$ for an
arbitrary set function $v$.

######  Lemma 3.1.

[32] Given any set function $v$, it can be expressed as a DS functions
$v(X)=f(X)-g(X),\forall X\subseteq V$ for some submodular functions $f$ and
$g$.

###### Proof.

Given a set function $v$, we can define $\alpha=\min_{X\subset Y\subseteq
V\setminus j}v(j|X)-v(j|Y)$111We denote $j,X,Y:X\subset Y\subseteq
V\setminus\\{j\\}$ by $X\subset Y\subseteq V\setminus j$.. Clearly $\alpha<0$,
since otherwise $v$ would be submodular. Now consider any (strictly)
submodular function $g$, i.e., one having $\beta=\min_{X\subset Y\subseteq
V\setminus j}g(j|X)-g(j|Y)>0$. Define
$f^{\prime}(X)=v(X)+\frac{|\alpha^{\prime}|}{\beta}g(X)$ with any
$\alpha^{\prime}\leq\alpha$. Now it is easy to see that $f^{\prime}$ is
submodular since $\min_{X\subset Y\subseteq V\setminus
j}f^{\prime}(j|X)-f^{\prime}(j|Y)\geq\alpha+|\alpha^{\prime}|\geq 0$. Hence
$v(X)=f^{\prime}(X)-\frac{|\alpha^{\prime}|}{\beta}g(X)$, is a difference
between two submodular functions. ∎

The above proof requires the computation of $\alpha$ and $\beta$ which has, in
general, exponential complexity. Using the construction above, however, it is
easy to find the decomposition $f$ and $g$ under certain conditions on $v$.

######  Lemma 3.2.

If $\alpha$ or at least a lower bound on $\alpha$ for any set function $v$ can
be computed in polynomial time, functions $f$ and $g$ corresponding to $v$ can
obtained in polynomial time.

###### Proof.

Define $g$ as $g(X)=\sqrt{|X|}$. Then $\beta=\min_{X\subset Y\subseteq
V\setminus j}\sqrt{|X|+1}-\sqrt{|X|}-\sqrt{|Y|+1}+\sqrt{|Y|}=\min_{X\subset
V\setminus
j}\sqrt{|X|+1}-\sqrt{|X|}-\sqrt{|X|+2}+\sqrt{|X|+1}=2\sqrt{n-1}-\sqrt{n}-\sqrt{n-2}$.
The last inequality follows since the smallest difference in gains will occur
at $|X|=n-2$. Hence $\beta$ is easily computed, and given a lower bound on
$\alpha$, from Lemma 3.1 the decomposition can be obtained in polynomial time.
A similar argument holds for $g$ being other concave functions over $|X|$. ∎

Algorithm 1 The submodular-supermodular (SubSup) procedure [32]

1: $X^{0}=\emptyset$ ; $t\leftarrow 0$ ;

2: while not converged (i.e., $(X^{t+1}\neq X^{t})$) do

3: Randomly choose a permutation $\sigma^{t}$ whose chain contains the set
$X^{t}$.

4: $X^{t+1}:=\operatorname*{argmin}_{X}f(X)-h^{g}_{X^{t},\sigma^{t}}(X)$

5: $t\leftarrow t+1$

6: end while

The submodular supermodular (SubSup) procedure is given in Algorithm 1. At
every step of the algorithm, we minimize a submodular function which can be
performed in strongly polynomial time [34, 37] although the best known
complexity is $O(n^{5}\eta+n^{6})$ where $\eta$ is the cost of a function
evaluation. Algorithm 1 is guaranteed to converge to a local minima and
moreover the algorithm monotonically decreases the function objective at every
iteration, as we show below.

######  Lemma 3.3.

[32] Algorithm 1 is guaranteed to decrease the objective function at every
iteration. Further, the algorithm is guaranteed to converge to a local minima
by checking at most $O(n)$ permutations at every iteration.

###### Proof.

The objective reduces at every iteration since:

$\displaystyle f(X^{t+1})-g(X^{t+1})$ $\displaystyle\overset{a}{\leq}$
$\displaystyle f(X^{t+1})-h^{g}_{X^{t},\sigma^{t}}(X^{t+1})$
$\displaystyle\overset{b}{\leq}$ $\displaystyle
f(X^{t})-h^{g}_{X^{t},\sigma^{t}}(X^{t})$ $\displaystyle\overset{c}{=}$
$\displaystyle f(X^{t})-g(X^{t})$

Where (a) follows since $h^{g}_{X^{t},\sigma^{t}}(X^{t+1})\leq g(X^{t+1})$,
and (b) follows since $X^{t+1}$ is the minimizer of
$f(X)-h^{g}_{X^{t},\sigma^{t}}(X)$, and (c) follows since
$h^{g}_{X^{t},\sigma^{t}}(X^{t})=g(X^{t})$ from the tightness of the modular
lower bound.

Further note that, if there is no improvement in the function value by
considering $O(n)$ permutations each with different elements at
$\sigma^{t}(|X^{t}|-1)$ and $\sigma^{t}(|X^{t}|+1)$, then this is equivalent
to a local minima condition on $v$ since
$h^{g}_{X^{t},\sigma^{t}}(S_{|X^{t}|+1}^{\sigma})=f(S_{|X^{t}|+1}^{\sigma})$
and
$h^{g}_{X^{t},\sigma^{t}}(S_{|X^{t}|-1}^{\sigma})=f(S_{|X^{t}|-1}^{\sigma})$.
∎

Algorithm 1 requires performing a submodular function minimization at every
iteration which while polynomial in $n$ is (due to the complexity described
above) not practical for large problem sizes. So while the algorithm reaches a
local minima, it can be costly to find it. A desirable result, therefore,
would be to develop new algorithms for minimizing DS functions, where the new
algorithms have the same properties as the SubSup procedure but are much
faster in practice. We give this in the following sections.

##  4 Alternate algorithms for minimizing DS functions

In this section we propose two new algorithms to minimize DS functions, both
of which are guaranteed to monotonically reduce the objective at every
iteration and converge to local minima. We briefly describe these algorithms
in the subsections below.

###  4.1 The supermodular-submodular (SupSub) procedure

In the submodular-supermodular procedure we iteratively minimized $f(X)-g(X)$
by replacing $g$ by it’s modular lower bound at every iteration. We can
instead replace $f$ by it’s modular upper bound as is done in Algorithm 2,
which leads to the supermodular-submodular procedure.

Algorithm 2 The supermodular-submodular (SupSub) procedure

1: $X^{0}=\emptyset$ ; $t\leftarrow 0$ ;

2: while not converged (i.e., $(X^{t+1}\neq X^{t})$) do

3: $X^{t+1}:=\operatorname*{argmin}_{X}m^{f}_{X^{t}}(X)-g(X)$

4: $t\leftarrow t+1$

5: end while

In the SupSub procedure, at every step we perform submodular maximization
which, although NP complete to solve exactly, admits a number of fast constant
factor approximation algorithms [4, 9]. Notice that we have two modular upper
bounds and hence there are a number of ways we can choose between them. One
way is to run both maximization procedures with the two modular upper bounds
at every iteration in parallel, and choose the one which is better. Here by
better we mean the one in which the function value is lesser. Alternatively we
can alternate between the two modular upper bounds by first maximizing the
expression using the first modular upper bound, and then maximize the
expression using the second modular upper bound. Notice that since we perform
approximate submodular maximization at every iteration, we are not guaranteed
to monotonically reduce the objective value at every iteration. If, however,
we ensure that at every iteration we take the next step only if the objective
$v$ does not increase, we will restore monotonicity at every iteration. Also,
in some cases we converge to local optima as shown in the following theorem.

######  Theorem 4.1.

Both variants of the supermodular-submodular procedure (Algorithm 2)
monotonically reduces the objective value at every iteration. Moreover,
assuming a submodular maximization procedure in line 3 that reaches a local
maxima of $m^{f}_{X^{t}}(X)-g(X)$, then if Algorithm 2 does not improve under
both modular upper bounds then it reaches a local optima of $v$.

###### Proof.

For either modular upper bound, we have:

$\displaystyle f(X^{t+1})-g(X^{t+1})$
$\displaystyle\overset{a}{\leq}m^{f}_{X^{t}}(X^{t+1})-g(X^{t+1})$
$\displaystyle\overset{b}{\leq}m^{f}_{X^{t}}(X^{t})-g(X^{t})$
$\displaystyle\overset{c}{=}f(X^{t})-g(X^{t}),$

where (a) follows since $f(X^{t+1})\leq m^{f}_{X^{t}}(X^{t+1})$, and (b)
follows since we assume that we take the next step only if the objective value
does not increase and (c) follows since $m^{f}_{X^{t}}(X^{t})=f(X^{t})$ from
the tightness of the modular upper bound.

To show that this algorithm converges to a local minima, we assume that the
submodular maximization procedure in line 3 converges to a local maxima. Then
observe that if the objective value does not decrease in an iteration under
both upper bounds, it implies that $m^{f}_{X^{t}}(X^{t})-g(X^{t})$ is already
a local optimum in that (for both upper bounds) we have
$m^{f}_{X^{t}}(X^{t}\cup j)-g(X^{t}\cup j)\geq
m^{f}_{X^{t}}(X^{t})-g(X^{t}),\forall j\notin X^{t}$ and
$m^{f}_{X^{t}}(X^{t}\backslash j)-g(X^{t}\backslash j)\geq
m^{f}_{X^{t}}(X^{t})-g(X^{t}),\forall j\in X^{t}$. Note that
$m^{f}_{X^{t},1}(X^{t}\backslash j)=f(X^{t})-f(j|X^{t}\backslash
j)=f(X^{t}\backslash j)$ and $m^{f}_{X^{t},2}(X^{t}\cup
j)=f(X^{t})+f(j|X^{t})=f(X^{t}\cup j)$ and hence if both modular upper bounds
are at a local optima, it implies
$f(X^{t})-g(X^{t})=m^{f}_{X^{t},1}(X^{t})-g(X^{t})\leq
m^{f}_{X^{t},1}(X^{t}\backslash j)-g(X^{t}\backslash j)=f(X^{t}\backslash
j)-g(X^{t}\backslash j)$. Similarly
$f(X^{t})-g(X^{t})=m^{f}_{X^{t},2}(X^{t})-g(X^{t})\leq
m^{f}_{X^{t},2}(X^{t}\cup j)-g(X^{t}\cup j)=f(X^{t}\cup j)-g(X^{t}\cup j)$.
Hence $X^{t}$ is a local optima for $v(X)=f(X)-g(X)$, since $v(X^{t})\leq
v(X^{t}\cup j)$ and $v(X^{t})\leq v(X^{t}\backslash j)$. ∎

To ensure that we take the largest step at each iteration, we can use the
recently proposed tight (1/2)-approximation algorithm in [4] for unconstrained
non-monotone submodular function maximization — this is the best possible in
polynomial time for the class of submodular functions independent of the P=NP
question. The algorithm is a form of bi-directional randomized greedy
procedure and, most importantly for practical considerations, is linear time
[4]. In practice we just use a combination of a form of a simple greedy
procedure, and the bi-directional randomized algorithm, by picking the best
amongst the two at every iteration. Since the randomized greedy algorithm is
$1/2$ approximate, the combination of the two procedures also will be $1/2$
approximate.

Lastly, note that this algorithm is closely related to a local search
heuristic for submodular maximization [9]. In particular, if instead of using
the greedy algorithm entirely at every iteration, we take only one local step,
we get a local search heuristic. Hence, via the SupSub procedure, we may take
larger steps at every iteration as compared to a local search heuristic.

###  4.2 The modular-modular (ModMod) procedure

The submodular-supermodular procedure and the supermodular-submodular
procedure were obtained by replacing $g$ by it’s modular lower bound and $f$
by it’s modular upper bound respectively. We can however replace both of them
by their respective modular bounds, as is done in Algorithm 3.

Algorithm 3 Modular-Modular (ModMod) procedure

1: $X^{0}=\emptyset$; $t\leftarrow 0$ ;

2: while not converged (i.e., $(X^{t+1}\neq X^{t})$) do

3: Choose a permutation $\sigma^{t}$ whose chain contains the set $X^{t}$.

4:
$X^{t+1}:=\operatorname*{argmin}_{X}m^{f}_{X^{t}}(X)-h^{g}_{X^{t},\sigma^{t}}(X)$

5: $t\leftarrow t+1$

6: end while

In this algorithm at every iteration we minimize only a modular function which
can be done in $O(n)$ time, so this is extremely easy (i.e., select all
negative elements for the smallest minimum, or all non-positive elements for
the largest minimum). Like before, since we have two modular upper bounds, we
can use any of the variants discussed in the subsection above. Moreover, we
are still guaranteed to monotonically decrease the objective at every
iteration and converge to a local minima.

######  Theorem 4.2.

Algorithm 3 monotonically decreases the function value at every iteration. If
the function value does not increase on checking $O(n)$ different permutations
with different elements at adjacent positions and with both modular upper
bounds, then we have reached a local minima of $v$.

###### Proof.

Again we can use similar reasoning as the earlier proofs and observe that:

$\displaystyle f(X^{t+1})-g(X^{t+1})$ $\displaystyle\leq
m^{f}_{X^{t}}(X^{t+1})-h^{g}_{X^{t},\sigma^{t}}(X^{t+1})$ $\displaystyle\leq
m^{f}_{X^{t}}(X^{t})-h^{g}_{X^{t},\sigma^{t}}(X^{t})$
$\displaystyle=f(X^{t})-g(X^{t})$

We see that considering $O(n)$ permutations each with different elements at
$\sigma^{t}(|X^{t}|-1)$ and $\sigma^{t}(|X^{t}|+1)$, we essentially consider
all choices of $g(X^{t}\cup j)$ and $g(X^{t}\backslash j)$, since
$h^{g}_{X^{t},\sigma^{t}}(S_{|X^{t}|+1})=f(S_{|X^{t}|+1})$ and
$h^{g}_{X^{t},\sigma^{t}}(S_{|X^{t}|-1})=f(S_{|X^{t}|-1})$. Since we consider
both modular upper bounds, we correspondingly consider every choice of
$f(X^{t}\cup j)$ and $f(X^{t}\backslash j)$. Note that at convergence we have
that $m^{f}_{X^{t}}(X^{t})-h^{g}_{X^{t},\sigma^{t}}(X^{t})\leq
m^{f}_{X^{t}}(X)-h^{g}_{X^{t},\sigma^{t}}(X),\forall X\subseteq V$ for $O(n)$
different permutations and both modular upper bounds. Correspondingly we are
guaranteed that (since the expression is modular) $\forall j\notin
X^{t},v(j|X^{t})\geq 0$ and $\forall j\in X^{t},v(j|X^{t}\backslash j)\geq 0$,
where $v(X)=f(X)-g(X)$. Hence the algorithm converges to a local minima. ∎

An important question is the choice of the permutation $\sigma^{t}$ at every
iteration $X^{t}$. We observe experimentally that the quality of the algorithm
depends strongly on the choice of permutation. Observe that $f(X)-g(X)\leq
m^{f}_{X^{t}}(X)-h^{g}_{X^{t},\sigma^{t}}(X)$, and
$f(X^{t})-g(X^{t})=m^{f}_{X^{t}}(X^{t})-h^{g}_{X^{t},\sigma^{t}}(X^{t})$.
Hence, we might obtain the greatest local reduction in the value of $v$ by
choosing permutation
$\sigma^{*}\in\operatorname*{argmin}_{\sigma}\min_{X}(m^{f}_{X^{t}}(X)-h^{g}_{X^{t},\sigma^{t}}(X))$,
or the one which maximizes $h^{g}_{X^{t},\sigma^{t}}(X)$. We in fact might
expect that choosing $\sigma^{t}$ ordered according to greatest gains of $g$,
with respect to $X^{t}$, we would achieve greater descent at every iteration.
Another choice is to choose the permutation $\sigma$ based on the ordering of
gains of $v$ (or even $m^{f}_{X^{t}}$). Through the former we are guaranteed
to at least progress as much as the local search heuristic. Indeed, we observe
in practice that the first two of these heuristics performs much better than a
random permutation for both the ModMod and the SubSup procedure, thus
addressing a question raised in [32] about which ordering to use. Practically
for the feature selection problem, the second heuristic seems to work the
best.

###  4.3 Constrained minimization of a difference between submodular functions

In this section we consider the problem of minimizing the difference between
submodular functions subject to constraints. We first note that the problem of
minimizing a submodular function under even simple cardinality constraints in
NP hard and also hard to approximate [38]. Since there does not yet seem to be
a reasonable algorithm for constrained submodular minimization at every
iteration, it is unclear how we would use Algorithm 1. However the problem of
submodular maximization under cardinality, matroid, and knapsack constraints
though NP hard admits a number of constant factor approximation algorithms
[33, 28] and correspondingly the cardinality constraints can be easily
introduced in Algorithm 2. Moreover, since a non-negative modular function can
be easily, directly and even exactly optimized under cardinality, knapsack and
matroid constraints [16], Algorithm 3 can also easily be utilized. In
addition, since problems such as finding the minimum weight spanning tree,
min-cut in a graph, etc., are polynomial time algorithms in a number of cases,
Algorithm 3 can be used when minimizing a non-negative function $v$
expressible as a difference between submodular functions under combinatorial
constraints. If $v$ is non-negative, then so is its modular upper bound, and
then the ModMod procedure can directly be used for this problem — each
iteration minimizes a non-negative modular function subject to combinatorial
constraints which is easy in many cases [16, 15].

##  5 Theoretical results

In this section we analyze the computational and approximation bounds for this
problem. For simplicity we assume that the function $v$ is normalized, i.e
$v(\emptyset)=0$. Hence we assume that $v$ achieves it minima at a negative
value and correspondingly the approximation factor in this case will be less
than $1$.

We note in passing that the results in this section are mostly negative, in
that they demonstrate theoretically how complex a general problem such as
$\min_{X}[f(X)-g(X)]$ is, even for submodular $f$ and $g$. In this paper,
rather than consider these hardness results pessimistically, we think of them
as providing justification for the heuristic procedures given in Section 4 and
[32]. In many cases, inspired heuristics can yield good quality and hence
practically useful algorithms for real-world problems. For example, the ModMod
procedure (Algorithm 3) and even the SupSub procedure (Algorithm 2) can scale
to very large problem sizes, and thus can provide useful new strategies for
the applications listed in Section 1.

###  5.1 Hardness

Observe that the class of DS functions is essentially the class of general set
functions, and hence the problem of finding optimal solutions is NP-hard. This
is not surprising since general set function minimization is inapproximable
and there exist a large class of functions where all (adaptive, possibly
randomized) algorithms perform arbitrarily poorly in polynomial time [39].
Clearly as is evident from Theorem 3.1, even the problem of finding the
submodular functions $f$ and $g$ requires exponential complexity. We moreover
show in the following theorem, however, that this problem is multiplicatively
inapproximable even when the functions $f$ and $g$ are easy to find.

######  Theorem 5.1.

Unless P = NP, there cannot exist any polynomial time approximation algorithm
for $\min_{X}v(X)$ where $v(X)=[f(X)-g(X)]$ is a positive set function and $f$
and $g$ are given submodular functions. In particular, let $n$ be the size of
the problem instance, and $\alpha(n)>0$ be any positive polynomial time
computable function of $n$. If there exists a polynomial-time algorithm which
is guaranteed to find a set
$X^{\prime}:f(X^{\prime})-g(X^{\prime})<\alpha(n)\mbox{OPT}$, where
OPT=$\min_{X}f(X)-g(X)$, then P = NP.

###### Proof.

We prove this by reducing this to the subset sum problem. Given a positive
modular function $m$ and a positive constant $t$, is there a subset
$S\subseteq V$ such that $m(S)=t$? First we choose a random set $C$ (unknown
to the algorithm), and define $t=m(C)$. Define a set function $v$, such that
$v(S)=1$, if $m(S)=t$ and $v(S)=\frac{1}{\alpha(n)}-o(1)$ otherwise. Observe
that $\min_{S}v(S)=\frac{1}{\alpha(n)}-o(1)$, since $\alpha(n)>1$. Note that
$\alpha=\min_{X\subset Y\subseteq V\setminus j}v(j|X)-v(j|Y)\geq
2(\frac{1}{\alpha(n)}-1)$. Hence we can easily compute a lower bound on
$\alpha$ and hence from lemma 3.2 we can directly compute the decomposition
$f$ and $g$. In fact notice that the decomposition is directly computable
since both $\alpha$ and $\beta$ are known.

Now suppose there exists a polynomial time algorithm for this problem with an
approximation factor of $\alpha(n)$. This implies that the algorithm is
guaranteed to find a set $S$, such that $v(S)<1$. Hence this algorithm will
solve the subset sum problem in polynomial time, which is a contradiction
unless P = NP. ∎

In fact we show below that independent of the $P=NP$ question, there cannot
exist a sub-exponential time algorithm for this problem with any constant
factor approximation. The theorem below gives information theoretic hardness
for this problem.

######  Theorem 5.2.

For any $0<\epsilon<1$, there cannot exist any deterministic (or possibly
randomized) algorithm for $\min_{X}[f(X)-g(X)]$ (where $f$ and $g$ are given
submodular functions), that always finds a solution which is at most
$\frac{1}{\epsilon}$ times the optimal, in fewer than $e^{\epsilon^{2}n/8}$
queries.

###### Proof.

For showing this theorem, we use the same proof technique as in [9]. Define
two sets $C$ and $D$, such that $V=C\cup D$ and $|C|=|D|=n/2$. We then define
a set function $v(S)$ which depends only on $k=|S\cap C|$ and $l=|S\cap D|$.
In particular define $v(S)=\frac{1}{\epsilon},\mbox{ if }|k-l|\leq\epsilon n$
and $v(S)=1,\mbox{ if }|k-l|>\epsilon n$. Again, we have a trivial bound on
$\alpha$ here since $v(j|X)\geq\frac{1}{\epsilon}-1$ and $v(j|Y)\leq
1-\frac{1}{\epsilon}$. Hence, $\alpha=\min_{X\subset Y\subseteq V\setminus
j}v(j|X)-v(j|Y)>2|1-\frac{1}{\epsilon}|$. Thus, for this set function, a
decomposition $v=f-g$ can easily be obtained (Lemma 3.2).

Now, let the partition $(C,D)$ be taken uniformly at random and unknown to the
algorithm. The algorithm issues some queries $S$ to the value oracle. Call $S$
“unbalanced” if $|S\cap C|$ differs from $|S\cap D|$ by more than $\epsilon
n$. Recall the Chernoff bounds [1]: Let $Y_{1},Y_{2},\cdots,Y_{t}$ be
independent random variables in $[-1,1]$, such that $\mathbb{E}[Y_{i}]=0$,
then:

$Pr[\sum_{i=1}^{t}Y_{i}>\lambda]\leq 2e^{-\lambda^{2}/2t}.$ (6)

Define $Y_{i}=I(i\in S)[I(i\in C)-I(i\in D)]$. Clearly $Y_{i}\in[-1,1]$, and
we can use the bounds above. Hence for any query $S$, the probability that $S$
is unbalanced is at most $2e^{-\epsilon^{2}n/2}$. Thus, we can see that even
after $e^{\epsilon^{2}n/4}$ number of queries, the probability that the
resulting set is unbalanced is still $2e^{-\epsilon^{2}n/4}$. Hence any
algorithm will query only balanced sets regardless of $C$ and $D$, and
consequently with high probability the algorithm will obtain
$\frac{1}{\epsilon}$ as the minimum, while the actual minimum is $1$. Thus,
such an algorithm will never be able to achieve an approximation factor better
than $\frac{1}{\epsilon}$. ∎

Essentially the theorems above say that even when we are given (or can easily
find) a decomposition such that $v(X)=f(X)-g(X)$, there exist set functions
such that any algorithm (either adaptive or randomized) cannot be approximable
upto any constant factor. It is possible that one could come up with an
information theoretic construction to show this same result for any polynomial
approximation factor. However under the assumption of P$\neq$NP, Theorem 5.1
shows that this problem is inapproximable upto any polynomial factor. Hence
any algorithm trying to find the global optimum for this problem [5, 20] can
only be exponential in the worst case.

Interestingly, the hardness results above holds even when the submodular
functions $f$ and $g$ are monotone. This follows from the following Lemma:

######  Lemma 5.1.

Given (not necessarily monotone) submodular functions $f$ and $g$, there
exists monotone submodular functions $f^{\prime}$ and $g^{\prime}$ such that,

$\displaystyle f(X)-g(X)=f^{\prime}(X)-g^{\prime}(X),\forall X\subseteq V$ (7)

###### Proof.

The proof of this Lemma follows from a simple observation. The decomposition
theorem of [7] shows that any submodular function can be decomposed into a
modular function plus a monotone non-decreasing and totally normalized
polymatroid rank function. Specifically, given submodular $f,g$ we have

$\displaystyle f^{\prime}(X)\triangleq f(X)-\sum_{j\in X}f(j|V\backslash j)$
(8)

and

$\displaystyle g^{\prime}(X)\triangleq g(X)-\sum_{j\in X}g(j|V\backslash j)$
(9)

$f^{\prime},g^{\prime}$ are then totally normalized polymatroid rank
functions. Hence we have: $v(X)=f^{\prime}(X)-g^{\prime}(X)+k(X)$, with
modular $k(X)=\sum_{j\in X}v(j|V\backslash j)$. The idea is then to add $v(j)$
to $f^{\prime}$ if $v(j)\geq 0$ or add it to $g^{\prime}$ other-wise. In
particular, let $V^{+}=\\{j:v(j)\geq 0\\}$ and $V^{-}=\\{j:v(j)<0\\}$. Notice
that $V^{+}\cup V^{-}=V$. Then,

$\displaystyle v(X)=f^{\prime}(X)+k(X\cap V^{+})-\\{g^{\prime}-k(X\cap
V^{-})\\}$ (10)

Notice above that $f^{\prime}(X)+k(X\cap V^{+})$ and $g^{\prime}-k(X\cap
V^{-})$ are both monotone non-decreasing. Hence proved. ∎

This then implies the following corollary.

######  Corollary 5.3.

Given submodular functions $f$ and $g$ such that $v(X)=f(X)-g(X)\geq 0$, the
problem $\min_{X\subseteq V}v(X)$ is inapproximable, even if both $f$ and $g$
are monotone non-decreasing submodular.

###  5.2 Polynomial time lower and upper bounds

Since any submodular function can be decomposed into a modular function plus a
monotone non-decreasing and totally normalized polymatroid rank function [7],
we have: $v(X)=f^{\prime}(X)-g^{\prime}(X)+k(X)$, with modular
$k(X)=\sum_{j\in X}v(j|V\backslash j)$ and $f^{\prime}$ and $g^{\prime}$ being
the totally normalized polymatroid functions.

The algorithms in the previous sections are all based on repeatedly finding
upper bounds for $v$. The following lower bounds directly follow from the
results above.

######  Theorem 5.4.

We have the following two lower bounds on the minimizers of $v(X)=f(X)-g(X)$:

$\displaystyle\min_{X}v(X)$
$\displaystyle\geq\min_{X}f^{\prime}(X)+k(X)-g^{\prime}(V)$
$\displaystyle\min_{X}v(X)$ $\displaystyle\geq
f^{\prime}(\emptyset)-g^{\prime}(V)+\sum_{j\in V}\min(k(j),0)$

###### Proof.

Notice that

$\displaystyle\min_{X}f(X)-g(X)$
$\displaystyle=\min_{X}f^{\prime}(X)-g^{\prime}(X)+k(X)$
$\displaystyle\geq\min_{X}(f^{\prime}(X)+k(X))-\max_{X}g^{\prime}(X)$
$\displaystyle=\min_{X}f^{\prime}(X)+k(X)-g^{\prime}(V)$

To get the second result, we start from the bound above and loosen it as:

$\displaystyle\min_{X}f^{\prime}(X)+k(X)-g^{\prime}(V)$
$\displaystyle\geq\min_{X}f^{\prime}(X)+\min_{X}k(X)-g^{\prime}(V)$
$\displaystyle=f^{\prime}(\emptyset)+\sum_{j\in V}\min(v(j|V\backslash
j),0)-g^{\prime}(V)$ $\displaystyle=f^{\prime}(\emptyset)+\sum_{j\in
V}\min(k(j),0)-g^{\prime}(V)$ (11)

∎

The above lower bounds essentially provide bounds on the minima of the
objective and thus can be used to obtain an additive approximation guarantee.
The algorithms described in this paper are all polynomial time algorithms (as
we show below) and correspondingly from the bounds above we can get an
estimate on how far we are from the optimal.

###  5.3 Computational Bounds

We now provide computational bounds for $\epsilon$-approximate versions of our
algorithms. Note that this was left as an open question in [32]. Finding the
local minimizer of DS functions is PLS complete since it generalizes the
problem of finding the local optimum of the MAX-CUT problem [36]. Note that
this trivially generalizes the MAX-CUT problem since if we set $f(X)=0$ and
$g(X)$ is the cut function, we get the max cut problem. However we show that
an $\epsilon$-approximate version of this algorithm will converge in
polynomial time.

######  Definition 5.1.

An $\epsilon$-approximate version of an iterative monotone non-decreasing
algorithm for minimizing a set function $v$ is defined as a version of that
algorithm, where we proceed to step $t+1$ only if $v(X^{t+1})\leq
v(X^{t})(1+\epsilon)$.

Note that the $\epsilon$-approximate versions of algorithms 1, 2 and 3, are
guaranteed to converge to $\epsilon$-approximate local optima. An
$\epsilon$-approximate local optima of a function $v$ is a set $X$, such that
$v(X\cup j)\geq v(X)(1+\epsilon)$ and $v(X\backslash j)\geq v(X)(1+\epsilon)$.
W.l.o.g., assume that $X^{0}=\emptyset$. Then we have the following
computational bounds:

######  Theorem 5.5.

The $\epsilon$-approximate versions of algorithms 1, 2 and 3 have a worst case
complexity of $O(\frac{\log(|M|/|m|)}{\epsilon}T))$, where
$M=f^{\prime}(\emptyset)+\sum_{j\in V}\min(v(j|V\backslash
j),0)-g^{\prime}(V)$, $m=v(X^{1})$ and $O(T)$ is the complexity of every
iteration of the algorithm (which corresponds to respectively the submodular
minimization, maximization, or modular minimization in algorithms 1, 2 and
3)..

###### Proof.

Observe that $m=v(X^{1})\leq v(X^{0})=0$. Correspondingly if $v(X^{1})=0$, it
implies that the algorithm has converged, and cannot improve (since we are
assuming our algorithms are $\epsilon-$approximate. Hence in this case the
algorithm will converge in one iteration. Consider then the case of $m<0$.
Note also from Theorem 5.4 that $M=f^{\prime}(\emptyset)+\sum_{j\in
V}\min(v(j|V\backslash j),0)-g^{\prime}(V)<0$ and that $\min_{X}f(X)-g(X)\geq
M$. Since we are guaranteed to improve by a factor by at least $1+\epsilon$ at
every iteration we have that in $k$ iterations:
$|m|(1+\epsilon)^{k}\leq|M|\Rightarrow k=O(\frac{\log(|M|/|m|)}{\epsilon})$.
Also since we assume that the complexity at every iteration is $O(T)$ we get
the above result. ∎

Observe that for the algorithms we use, $O(T)$ is strongly polynomial in $n$.
The best strongly polynomial time algorithm for submodular function
minimization is $O(n^{5}\eta+n^{6})$ [34] (the lower bound is currently
unknown). Further the worst case complexity of the greedy algorithm for
maximization is $O(n^{2})$ while the complexity of modular minimization is
just $O(n)$. Note finally that these are worst case complexities and actually
the algorithms run much faster in practice.

(a) SVM

(b) NB

Figure 1: Plot showing the accuracy rates vs. the number of features on the
Mushroom data set.

##  6 Experiments

We test our algorithms on the feature subset selection problem in the
supervised setting. Given a set of features
$X_{V}=\\{X_{1},X_{2},\cdots,X_{|V|}\\}$, we try to find a subset of these
features $A$ which has the most information from the original set $X_{V}$
about a class variable $C$ under constraints on the size or cost of $A$.
Normally the number of features $|V|$ is quite large and thus the training and
testing time depend on $|V|$. In many cases, however, there is a strong
correlation amongst features and not every feature is novel. We can thus
perform training and testing with a much smaller number of features $|A|$
while obtaining (almost) the same error rates.

The question is how to find the most representative set of features $A$. The
mutual information between the chosen set of features and the target class
$C$, $I(X_{A};C)$, captures the relevance of the chosen subset of features. In
most cases the selected features are not independent given the class $C$ so
the naïve Bayes assumption is not applicable, meaning this is not a pure
submodular optimization problem. As mentioned in Section 1, $I(X_{A};C)$ can
be exactly expressed as a difference between submodular functions $H(X_{A})$
and $H(X_{A}|C)$.

###  6.1 Modular Cost Feature Selection

In this subsection, we look at the problem of maximizing
$I(X_{A};C)-\lambda|A|$, as a regularized feature subset selection problem.
Note that a mutual information $I(X_{A};C)$ query can easily be estimated from
the data by just a single sweep through this data. Further we have observed
that using techniques such as Laplace smoothing helps to improve mutual
information estimates without increasing computation. In these experiments,
therefore, we estimate the mutual information directly from the data and run
our algorithms to find the representative subset of features.

We compare our algorithms on two data sets, i.e., the Mushroom data set [14]
and the Adult data set [22] obtained from [10]. The Mushroom data set has 8124
examples with 112 features, while the Adult data set has 32,561 examples with
123 features. In our experiments we considered subsets of features of sizes
between 5%-20% of the total number of features by varying $\lambda$. We tested
the following algorithms for the feature subset selection problem. We
considered two formulations of the mutual information, one under naïve Bayes,
where the conditional entropy $H(X_{A}|C)$ can be written as
$H(X_{A}|C)=\sum_{j\in A}H(X_{i}|C)$ and another where we do not assume such
factorization. We call these two formulations factored and non-factored
respectively. We then considered the simple greedy algorithm, of iteratively
adding features at every step to the factored and non-factored mutual
information, which we call GrF and GrNF respectively. Lastly, we use the new
algorithms presented in this paper on the non-factored mutual information.

We then compare the results of the greedy algorithms with those of the three
algorithms for this problem, using two pattern classifiers based on either a
linear kernel SVM (using [6]) or a naïve Bayes (NB) classifier. We call the
results obtained from the supermodular-submodular heuristic as “SupSub”, the
submodular-supermodular procedure [32] as “SubSup”, and the modular-modular
objective as “ModMod.” In the SubSup procedure, we use the minimum norm point
algorithm [12] for submodular minimization, and in the SubSup procedure, we
use the optimal algorithm of [4] for submodular maximization. We observed that
the three heuristics generally outperformed the two greedy procedures, and
also that GRF can perform quite poorly, thus justifying our claim that the
naïve Bayes assumption can be quite poor. This also shows that although the
greedy algorithm in that case is optimal, the features are correlated given
the class and hence modeling it as a difference between submodular functions
gives the best results. We also observed that the SupSub and ModMod procedures
perform comparably to the SubSup procedure, while the SubSup procedure is
_much_ slower in practice. Comparing the running times, the ModMod and the
SupSub procedure are each a few times slower then the greedy algorithm (ModMod
is slower due computing the modular semigradients), while the SubSup procedure
is around 100 times slower. The SubSup procedure is slower due to general
submodular function minimization which can be quite slow.

The results for the Mushroom data set are shown in Figure 1. We performed a 10
fold cross-validation on the entire data set and observed that when using all
the features SVM gave an accuracy rate of 99.6% while the all-feature NB model
had an accuracy rate of 95.5%. The results for the Adult database are in
Figure 2. In this case with the entire set of features the accuracy rate of
SVM on this data set is 83.9% and NB is 82.3%.

(a) SVM

(b) NB

Figure 2: Plot showing the accuracy rates vs. the number of features on the
Adult data set.

In the mushroom data, the SVM classifier significantly outperforms the NB
classifier and correspondingly GrF performs much worse than the other
algorithms. Also, in most cases the three algorithms outperform GrNF. In the
adult data set, both the SVM and NB perform comparably although SVM
outperforms NB. However in this case also we observe that our algorithms
generally outperform GrF and GrNF.

###  6.2 Submodular cost feature selection

(a) SVM

(b) NB

Figure 3: Plot showing the accuracy rates vs. the cost of features for the
Mushroom data set

(a) SVM

(b) NB

Figure 4: Plot showing the accuracy rates vs. the cost of features for the
Adult data set

We perform synthetic experiments for the feature subset selection problem
under submodular costs. The cost model we consider is
$c(A)=\sum_{i}\sqrt{m(A\cap S_{i})}$. We partitioned $V$ into sets
$\\{S_{i}\\}_{i}$ and chose the modular function $m$ randomly. In this set of
experiments, we compare the accuracy of the classifiers vs. the cost
associated with the choice of features for the algorithms. Recall, with simple
(modular) cardinality costs the greedy algorithms performed decently in
comparison to our algorithms in the adult data set, where the NB assumption is
reasonable. However with submodular costs, the objective is no longer
submodular even under the NB assumption and thus the greedy algorithms perform
much worse. This is unsurprising since the greedy algorithm is approximately
optimal only for monotone submodular functions. This is even more strongly
evident from the results of the mushrooms data-set (Figure 3)

##  7 Discussion

We have introduced new algorithms for optimizing the difference between two
submodular functions, provided new theoretical understanding that provides
some justification for heuristics, have outlined applications that can make
use of our procedures, and have tested in the case of feature selection with
modular and submodular cost features. Our new ModMod procedure is fast at each
iteration and experimentally does about as well as the SupSub and SubSup
procedures. The ModMod procedure, moreover, can also be used under various
combinatorial constraints, and therefore the ModMod procedure may hold the
greatest promise as a practical heuristic. An alternative approach, not yet
evaluated, would be to try the convex-concave procedure [41] on the Lovász
extensions of $f$ and $g$ since subgradients in such case are so easy to
obtain.

Acknowledgments: We thank Andrew Guillory, Manas Joglekar, Stefanie Jegelka,
and the rest of the submodular group at UW for discussions. This material is
based upon work supported by the National Science Foundation under Grant No.
(IIS-1162606), and is also supported by a Google, a Microsoft, and an Intel
research award.

## References

  * Alon and Spencer [2000] N. Alon and J. H. Spencer.  The probabilistic method, 2nd ed.  _Wiley-Interscience, New York_ , 2000.
  * Boros and Hammer [2002] E. Boros and P. L. Hammer.  Pseudo-boolean optimization.  _Discrete Applied Math._ , 123(1–3):155 – 225, 2002.  ISSN 0166-218X.  doi: 10.1016/S0166-218X(01)00341-9.  URL http://www.sciencedirect.com/science/article/pii/S0166218X01003419.
  * Brenner et al. [2000] N. Brenner, S.P. Strong, R. Koberle, W. Bialek, and R.R.R. Steveninck.  Synergy in a neural code.  _Neural Computation_ , 12(7):1531–1552, 2000\.
  * Buchbinder et al. [2012] N. Buchbinder, M. Feldman, J. Naor, and R. Schwartz.  A tight (1/2) linear-time approximation to unconstrained submodular maximization.  _In FOCS_ , 2012.
  * Byrnes [2009] K. Byrnes.  Maximizing general set functions by submodular decomposition.  _Arxiv preprint arXiv:0906.0120_ , 2009.
  * Chang and Lin [2011] Chih-Chung Chang and Chih-Jen Lin.  LIBSVM: A library for support vector machines.  _ACM Transactions on Intelligent Systems and Technology_ , 2, 2011\.
  * Cunningham [1983] W. H. Cunningham.  Decomposition of submodular functions.  _Combinatorica_ , 3(1):53–68, 1983.
  * Edmonds [1970] J. Edmonds.  Submodular functions, matroids and certain polyhedra.  _Combinatorial structures and their Applications_ , 1970.
  * Feige et al. [2011] Uriel Feige, Vahab Mirrokni, and Jan Vondrák.  Maximizing non-monotone submodular functions.  _SIAM J. COMPUT._ , 40(4):1133–1155, 2011.
  * Frank and Asuncion [2010] A. Frank and A. Asuncion.  UCI machine learning repository, 2010.  URL http://archive.ics.uci.edu/ml.
  * Fujishige [2005] S. Fujishige.  _Submodular functions and optimization_ , volume 58.  Elsevier Science, 2005.
  * Fujishige and Isotani [2011] S. Fujishige and S. Isotani.  A submodular function minimization algorithm based on the minimum-norm base.  _Pacific Journal of Optimization_ , 7:3–17, 2011.
  * Gallagher et al. [2011] Andrew C Gallagher, Dhruv Batra, and Devi Parikh.  Inference for order reduction in markov random fields.  In _Computer Vision and Pattern Recognition (CVPR), 2011 IEEE Conference on_ , pages 1857–1864. IEEE, 2011.
  * Iba et al. [1988] W. Iba, J. Wogulis, and P. Langley.  Trading off simplicity and coverage in incremental concept learning.  In _Proceedings of Fifth International Conference on Machine Learning_ , pages 73–79, 1988.
  * Jegelka and Bilmes [2010] S. Jegelka and J. Bilmes.  Cooperative cuts: Graph cuts with submodular edge weights.  Technical report, Technical Report TR-189, Max Planck Institute for Biological Cybernetics, 2010.
  * Jegelka and Bilmes [2011a] S. Jegelka and J. Bilmes.  Online submodular minimization for combinatorial structures.  _ICML_ , 2011a.
  * Jegelka and Bilmes [2011b] S. Jegelka and J. Bilmes.  Submodularity beyond submodular energies: coupling edges in graph cuts.  In _Computer Vision and Pattern Recognition (CVPR)_ , 2011b.
  * Jegelka and Bilmes [2011c] S. Jegelka and J. A. Bilmes.  Approximation bounds for inference using cooperative cuts.  In _ICML_ , 2011c.
  * Jegelka and Bilmes [2011d] S. Jegelka and J. A. Bilmes.  Submodularity beyond submodular energies: coupling edges in graph cuts.  In _CVPR_ , 2011d.
  * Kawahara and Washio [2011] Yoshinobu Kawahara and Takashi Washio.  Prismatic algorithm for discrete dc programming problems.  In _NIPS_ , 2011.
  * Kempe et al. [2003] D. Kempe, J. Kleinberg, and E. Tardos.  Maximizing the spread of influence through a social network.  In _9th ACM SIGKDD Intl. Conf. on Knowledge Discovery and Data Mining_ , 2003.
  * Kohavi [1996] R. Kohavi.  Scaling up the accuracy of naive-bayes classifiers: A decision-tree hybrid.  In _Proceedings of the second international conference on knowledge discovery and data mining_ , volume 7, 1996.
  * Kolmogorov and Rother [2007] V. Kolmogorov and C. Rother.  Minimizing nonsubmodular functions with graph cuts–a review.  _IEEE TPAMI_ , 29(7):1274–1279, 2007.
  * Kolmogorov and Zabih [2004] V. Kolmogorov and R. Zabih.  What energy functions can be minimized via graph cuts?  _IEEE TPAMI_ , 26(2):147–159, 2004.
  * Krause and Guestrin [2005] A. Krause and C. Guestrin.  Near-optimal nonmyopic value of information in graphical models.  In _Proceedings of Uncertainity in Artificial Intelligence_. UAI, 2005.
  * Krause et al. [2008a] A. Krause, A. Singh, and C. Guestrin.  Near-optimal sensor placements in Gaussian processes: Theory, efficient algorithms and empirical studies.  _JMLR_ , 9:235–284, 2008a.
  * Krause et al. [2008b] Andreas Krause, Brendan McMahan, Carlos Guestrin, and Anupam Gupta.  Robust submodular observation selection.  _Journal of Machine Learning Research (JMLR)_ , 9:2761–2801, 2008b.
  * Lee et al. [2009] J. Lee, V.S. Mirrokni, V. Nagarajan, and M. Sviridenko.  Non-monotone submodular maximization under matroid and knapsack constraints.  In _STOC_ , pages 323–332. ACM, 2009.
  * Lin and Bilmes [2010] H. Lin and J. Bilmes.  Multi-document summarization via budgeted maximization of submodular functions.  _In NAACL_ , 2010.
  * Lin and Bilmes [2011a] H. Lin and J. Bilmes.  A class of submodular functions for document summarization.  _In ACL_ , 2011a.
  * Lin and Bilmes [2011b] H. Lin and J. A. Bilmes.  Optimal selection of limited vocabulary speech corpora.  In _Annual Conference of the International Speech Communication Association (INTERSPEECH)_ , Florence, Italy, 2011b.
  * Narasimhan and Bilmes [2005] M. Narasimhan and J. Bilmes.  A submodular-supermodular procedure with applications to discriminative structure learning.  In _UAI_ , 2005.
  * Nemhauser et al. [1978] G.L. Nemhauser, L.A. Wolsey, and M.L. Fisher.  An analysis of approximations for maximizing submodular set functions—i.  _Mathematical Programming_ , 14(1):265–294, 1978\.
  * Orlin [2009] J.B. Orlin.  A faster strongly polynomial time algorithm for submodular function minimization.  _Mathematical Programming_ , 118(2):237–251, 2009\.
  * Pearl [1988] J. Pearl.  _Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference_.  Morgan Kaufmann, 2nd printing edition, 1988.
  * Schäffer [1991] A.A. Schäffer.  Simple local search problems that are hard to solve.  _SIAM journal on Computing_ , 20:56, 1991.
  * Schrijver [2000] A. Schrijver.  A combinatorial algorithm minimizing submodular functions in strongly polynomial time.  _Journal of Combinatorial Theory, Series B_ , 80(2):346–355, 2000.
  * Svitkina and Fleischer [2008] Z. Svitkina and L. Fleischer.  Submodular approximation: Sampling-based algorithms and lower bounds.  In _FOCS_ , pages 697–706, 2008.
  * Trevisan [2004] L. Trevisan.  Inapproximability of combinatorial optimization problems.  _The Computing Research Repository_ , 2004.
  * Wainwright and Jordan [2008] M.J. Wainwright and M.I. Jordan.  Graphical models, exponential families, and variational inference.  _Foundations and Trends® in Machine Learning_ , 1(1-2):1–305, 2008.
  * Yuille and Rangarajan [2002] A.L. Yuille and A. Rangarajan.  The concave-convex procedure (CCCP).  _In NIPS_ , 2002.