text:
# A proof of Parisi’s conjecture on the random assignment problem

## 1. Introduction

The problem of minimizing the sum of $`k`$ elements in a matrix of nonnegative real numbers under the condition that no two of them may be in the same row or column is called an *assignment problem*. A set of matrix positions no two in the same row or column is called an *independent* set. An independent set of $`k`$ matrix positions will also be called a $`k`$-assignment.

A *random assignment problem*, or RAP for short, is given by a number $`k`$, and an $`m`$ by $`n`$ matrix ($`\mathrm{min}(m,n)k`$) of random variables. If $`P`$ is a random assignment problem, we denote by $`E(P)`$ the expected value of the minimal sum of an independent set of $`k`$ matrix elements.

In this article we prove the following.

###### Theorem 1.1 (Parisi’s Conjecture \[P98\]).

Let $`P`$ be the RAP where $`k=m=n`$ and the matrix entries are independent exponential random variables with intensity 1. Then

$$E(P)=1+\frac{1}{4}+\frac{1}{9}+\mathrm{}+\frac{1}{k^2}.$$

We also prove the following two generalizations.

###### Theorem 1.2 (Conjectured by D. Coppersmith and G. B. Sorkin \[CS98\]).

Let $`P`$ be an RAP where the matrix entries are independent exponential random variables with intensity 1. Then

(1)
$$E(P)=\underset{i+j<k}{}\frac{1}{(mi)(nj)}.$$

###### Theorem 1.3 (Conjectured in \[LW00\]).

Let $`P`$ be an RAP where some matrix entries are zero and all the other entries are independent exponential random variables with intensity 1. Then

(2)
$$E(P)=\frac{1}{mn}\underset{i,j}{}\frac{d_{i,j}(P)}{\left(\genfrac{}{}{0pt}{}{m1}{i}\right)\left(\genfrac{}{}{0pt}{}{n1}{j}\right)},$$

where $`d_{i,j}(P)`$ is an integer coefficient defined in terms of the combinatorics of the set of zeros, see Section 2.2.

The identity (2) will be referred to as the *cover formula*. When $`P`$ has no zeros,

$$d_{i,j}(P)=\left(\genfrac{}{}{0pt}{}{m}{i}\right)\left(\genfrac{}{}{0pt}{}{n}{j}\right).$$

Hence Theorem 1.3 $``$ Theorem 1.2. To see the implication Theorem 1.2 $``$ Theorem 1.1 (proved in \[CS98\]), note that if we put $`m=n=k`$ in (1), then the terms for which $`\mathrm{gcd}(ki,kj)=d`$ sum to $`1/d^2`$ for $`d=1,\mathrm{},k`$.

As a consequence of Theorem 1.1 we also obtain a new and completely different proof of the following theorem, conjectured by M. Mézard and G. Parisi \[MP85\].

###### Theorem 1.4 (D. Aldous \[A92, A01\]).

Let $`P`$ be the RAP where $`k=m=n`$ and the matrix entries are independent exponential random variables with intensity 1. Then

$$\underset{k\mathrm{}}{lim}E(P)=\frac{\pi ^2}{6}.$$

In Section 7 we mention some other corollaries of Theorem 1.3.

### 1.1. Outline of the proof

A key result, proved in Section 4, is a formula for the probability that a nonzero row in an RAP is used in the optimal $`k`$-assignment, see Theorem 4.1. From \[LW00\] we know that the probability that a nonzero element in an RAP $`P`$ is used in the optimal $`k`$-assignment can be written $`E(P)E(P^{})`$, where $`P^{}`$ is obtained from $`P`$ by setting the matrix element in question to zero, see Theorem 2.2. Therefore, the formula for the probability that a row (or column) is used gives linear recursions for the values of the RAP’s. Provided $`m`$ or $`n`$ is sufficiently large compared to $`k`$, this system of linear recursions has a unique solution given by the cover formula in Theorem 1.3, see Section 5. Finally in Section 6 we prove that for fixed $`k`$, fixed $`m`$, and a fixed set of zeros, $`E(P)`$ is given by a rational function in $`n`$, which must then agree with the cover formula.

### 1.2. Background

Random assignment problems have attracted the attention of researchers from physics, optimization, and probability. There are experimental results in \[O92, PR93\]. Constructive upper and lower bounds on $`E(P)`$ have been given in \[W79, O92, CS98, L93, K87, GK93\]. Mézard and Parisi \[MP85\] used the non-rigorous *replica method* and arrived at the conjectured limit $`\pi ^2/6`$. This limit was subsequently established rigorously by Aldous \[A92, A01\] using the *weak convergence* method on a weighted infinite tree model. In this paper we continue the *exact formulas*-approach inspired by \[P98\] and developed further in \[AS02, BCR02, CS98, CS02, LW00, EES01\].

There are also interesting results on similar problems such as finding a minimal spanning tree in a graph with random edge weights \[BFM98, FM89, F85, EES01\]. An intriguing question is why the Riemann $`\zeta `$-function appears in the limit of both the spanning tree and bipartite matching problems.

One tool that we believe can be useful to a wider range of problems is Theorem 2.2 below and its generalization Theorem 7.3 of \[LW00\].

## 2. Preliminaries

### 2.1. Probabilistic preliminaries

We say that a random variable $`X`$ is exponentially distributed with intensity $`a`$ if $`Pr(X>t)=e^{at}`$ for $`t0`$. The intensity of $`X`$ is denoted $`I(X)`$. We have $`E(X)=1/I(X)`$.

A random assignment problem $`P`$ is called *standard* if the matrix entries are either zero or independent exponentially distributed with intensity 1. A standard RAP is determined by the numbers $`k`$, $`m`$ and $`n`$, and the set $`Z`$ of zero elements.

The following is a well-known lemma.

###### Lemma 2.1.

Let $`X_1,\mathrm{}X_n`$ be independent exponential variables with intensities $`a_1,\mathrm{},a_n`$ respectively. Then the probability that $`X_i`$ is minimal among $`X_1,\mathrm{}X_n`$ is

$$\frac{a_i}{a_1+\mathrm{}+a_n}.$$

The minimum is an exponential variable $`Y`$ of intensity $`a_1+\mathrm{}+a_n`$ which is independent of which variable is minimal. Under the condition that $`X_i`$ is minimal, $`X_i=Y`$, and for $`ji`$, we can write $`X_j=Y+X_j^{}`$, where $`X_j^{}`$ is exponential of intensity $`a_i`$, and the variables $`Y`$ and $`X_j`$ for $`ji`$ are all independent.

We say that a deterministic assignment problem is *generic* if no sum of a set of nonzero matrix elements is equal to the sum of a different set of nonzero elements. In the RAP’s that we consider, the distributions of the nonzero matrix elements are continuous. Hence all RAP’s considered here will be generic with probability 1. In the generic case, a nonzero element is used either in every optimal $`k`$-assignment, or in none. Hence we can without ambiguity speak of the probability that a certain nonzero element is used in *the* optimal $`k`$-assignment, without specifying whether we take this to mean some optimal $`k`$-assignment, or every optimal $`k`$-assignment.

Notice that this is not the case for zero elements. Even in a generic case, there may be several different optimal $`k`$-assignments, that differ in the choice of zero elements.

The following theorem, Theorem 2.10 of \[LW00\], is essential for the recursion equations in Section 5.

###### Theorem 2.2.

Suppose $`P`$ is a standard RAP. Suppose that the entry in position $`(i,j)`$ is not zero. Let $`P^{}`$ be the standard RAP where we have replaced the entry in position $`(i,j)`$ in $`P`$ with a zero. Then the probability that $`(i,j)`$ belongs to an optimal $`k`$-assignment in $`P`$ is $`E(P)E(P^{})`$.

### 2.2. Covers

We will consider sets of rows and columns in the $`m\times n`$-matrix. A set $`\lambda `$ of rows and columns is said to cover a set of zeros $`Z`$ if every matrix position in $`Z`$ is either in a row or in a column that belongs to $`\lambda `$. A cover with $`s`$ rows and columns will be called an $`s`$-cover. $`k1`$-covers will be of particular importance. By a partial $`k1`$-cover of $`Z`$, we mean a set of rows and columns which is a subset of a $`k1`$-cover of $`Z`$.

By a cover of an RAP we mean a cover of its set of zeros. The *cover coefficient* $`d_{i,j}(P)`$ is the number of partial $`k1`$-covers of $`P`$ with $`i`$ rows and $`j`$ columns. For this to be nonzero, $`i`$ and $`j`$ have to be nonnegative integers with $`i+j<k`$. It is convenient to regard the cover coefficient as well-defined, but zero, for integers $`i,j`$ outside this range.

We say that a set $`\lambda `$ of rows and columns is an *optimal cover* of $`P`$, if $`\lambda `$ covers $`P`$, and $`\lambda `$ has minimal cardinality among all covers of $`P`$. The following lemma is well-known. For a general introduction to matching theory we refer to \[LP86\].

###### Lemma 2.3 (Lattice property of optimal covers).

The set of optimal covers of a set of matrix positions forms a lattice, where one of the lattice operations is taking union of row sets and intersection of column sets, and the opposite lattice operation is taking intersection of row sets and union of column sets.

In particular, there is a row-maximal optimal cover containing every row that belongs to some optimal cover, and similarly a column-maximal optimal cover containing every column that belongs to some optimal cover.

An RAP can be reformulated in a setting of bipartite graphs with random weights on the edges. An assignment is then a matching. In this setting a cover of rows and columns is a vertex cover of the subgraph of edges with weight zero. We will borrow terminology from matching theory and speak of paths, etc. referring to the corresponding graph concepts.

The following two theorems were together with Theorem 2.2 our main tools for computing $`E(P)`$ recursively in \[LW00\].

###### Theorem 2.4.

Let $`A`$ be a deterministic assignment problem. If there is no set of $`k`$ independent zeros, then each row and column which belongs to an optimal cover of the zeros must be used in every optimal $`k`$-assignment.

Let $`c`$ be an optimal cover, and let $`x`$ be a positive real number smaller than or equal to the minimum of the elements not covered by $`c`$. Let $`A^{}`$ be obtained from $`A`$ by subtracting $`x`$ from the elements not covered by $`c`$, and adding $`x`$ to the doubly covered elements. Then $`x`$ has been subtracted from the optimal $`k`$-assignment $`k|c|`$ times, and consequently

$$E(A)=(k|c|)x+E(A^{}).$$

The previous theorem will be used for random assignment problems, by conditioning on where the minimal non-covered element is.

###### Theorem 2.5.

Suppose $`P`$ is a standard RAP. Suppose also that a column $`c`$ is used in every $`k1`$-cover of $`P`$. In particular, this is the case if at least $`k`$ elements in $`c`$ are zero. Let $`P^{}`$ be the RAP with column $`c`$ deleted from $`P`$ and $`k`$ and $`n`$ decreased by one. Then

(3)
$$E(P)=E(P^{}).$$

Similarly, if a row belongs to every $`k1`$-cover of $`P`$, it can be deleted.

## 3. Combinatorics of two optimal $`k`$-assignments

In this section we prove some results that are needed in Section 4. Consider a deterministic assignment problem. If there are two different $`k`$-assignments $`\mu `$ and $`\nu `$, their symmetric difference $`\mu \mathrm{}\nu :=(\mu \backslash \nu )(\nu \backslash \mu )`$ will be of special importance. By a *$`(\mu ,\nu )`$-alternating path* we mean a sequence $`(i_1,j_1),(i_2,j_2),\mathrm{},(i_r,j_r)`$ of matrix positions, where positions belonging to $`\mu `$ alternate with positions belonging to $`\nu `$ and $`i_x=i_{x+1}`$ or $`j_x=j_{x+1}`$, for all $`x=1,\mathrm{},r1`$. The parameter $`r`$ will be called the *length* of the path. It is easy to see that the positions in $`\mu \mathrm{}\nu `$ will form $`(\mu ,\nu )`$-alternating paths in the matrix.

If two $`k`$-assignments $`\mu `$ and $`\nu `$ differ only at zero positions we say that $`\mu `$ and $`\nu `$ are *equivalent*, $`\mu \nu `$.

A deterministic assignment problem is *semi-generic* if there is exactly one nontrivial sum of nonzero elements that equals a sum of a distinct set of nonzero elements.

###### Lemma 3.1.

Let $`A`$ be a deterministic assignment problem. Let $`\mu `$ be an optimal $`k`$-assignment, and let $`a`$ be an element which does not belong to $`\mu `$. Suppose that there is another optimal $`k`$-assignment $`\nu `$ that contains $`a`$. Then there is an optimal $`k`$-assignment $`\nu ^{}`$ containing $`a`$, such that $`\mu \mathrm{}\nu ^{}`$ consists either of one path, or of two paths of odd length.

###### Proof.

Let $`T`$ be the path of $`\mu \mathrm{}\nu `$ that contains $`a`$. If $`T`$ is of even length, then we let $`\nu ^{}=\mu \mathrm{}T`$. If $`T`$ has odd length, then it contains one more element of one of $`\mu `$ and $`\nu `$ than of the other. Since $`\mu `$ and $`\nu `$ have the same size, there must be another path $`T^{}`$ of $`\mu \mathrm{}\nu `$ that balances, so that $`TT^{}`$ has equally many elements from $`\mu `$ and $`\nu `$. Then we let $`\nu ^{}=\mu \mathrm{}(TT^{})`$. Clearly $`\nu ^{}`$ has the desired properties. ∎

###### Lemma 3.2.

Let $`A`$ be a semi-generic matrix with two inequivalent optimal $`k`$-assignments $`\mu `$ and $`\nu `$. Then $`\mu \mathrm{}\nu `$ has at most two paths with some nonzero element, and if there are two, they both have odd length.

###### Proof.

By Lemma 3.1 there is an optimal $`k`$-assignment $`\nu ^{}`$ that contains a nonzero element not in $`\mu `$, such that $`\mu \mathrm{}\nu ^{}`$ has either one path, or two paths of odd length. Since there are only two equivalence classes of optimal $`k`$-assignments, we must have $`\nu ^{}\nu `$. ∎

###### Lemma 3.3.

Let $`A`$ be a deterministic assignment problem. Let $`r_1`$ and $`r_2`$ be rows, and let $`c_1`$ and $`c_2`$ be columns. Suppose that $`A(r_1,c_1)=0`$, $`A(r_2,c_1)>0`$, and $`A(r_2,c_2)>0`$. Suppose further that there is an optimal $`k`$-assignment that uses $`(r_2,c_2)`$, but not row $`r_1`$. Then there is no optimal $`k`$-assignment that uses $`(r_2,c_1)`$.

###### Proof.

Let $`\mu `$ be an optimal $`k`$-assignment that uses $`(r_2,c_2)`$, and suppose that there is an optimal $`k`$-assignment $`\nu `$ that uses $`(r_2,c_1)`$. By Lemma 3.1, we may assume that $`\mu \mathrm{}\nu `$ consists of at most two paths, and that if there are two paths, both contain an odd number of matrix elements.

There must be a zero element $`(s_1,c_1)\mu `$, otherwise $`(r_2,c_2)`$ could be replaced by $`(r_1,c_1)`$ in $`\mu `$. Let $`S`$ be the path of $`(\mu \mathrm{}\nu )\{(r_2,c_1)\}`$ that contains $`(s_1,c_1)`$. In other words, $`S`$ is the set of matrix positions in the $`(\mu ,\nu )`$-alternating path that starts at $`(s_1,c_1)`$ and continues in the direction opposite to that of $`(r_2,c_1)`$.

There must be an element $`(r_1,d)`$ of $`\nu `$ since otherwise $`(r_2,c_1)`$ could be replaced by $`(r_1,c_1)`$ in $`\nu `$. Since $`\mu `$ does not use row $`r_1`$, the element $`(r_1,d)`$ is in a path of $`\mu \mathrm{}\nu `$ which is a path with one end belonging to $`\nu `$. If the other end belongs to $`\mu `$, the path has an even number of elements, and must therefore be the same as the path containing $`(r_2,c_1)`$. Otherwise the other end too belongs to $`\nu `$. Then the other path has both ends in $`\mu `$. In either case, $`S`$ must end with an element of $`\mu `$ in a row which is not used by $`\nu `$.

If the sum of the matrix entries in $`\mu S`$ is greater than the sum of the matrix entries in $`\nu S`$, then the cost of the $`k`$-assignment $`(\mu \mathrm{}S)\{(r_1,c_1)\}`$ is smaller than the cost of $`\mu `$, contradicting the optimality of $`\mu `$. Otherwise the sum of the matrix entries in $`\mu S`$ is smaller than or equal to the sum of the matrix entries in $`\nu S`$. Then the cost of $`(\nu \mathrm{}S)\{(r_2,c_1)\}`$ is smaller than the cost of $`\nu `$, again a contradiction. ∎

###### Lemma 3.4.

Let $`A`$ be a semi-generic deterministic assignment problem, and suppose that there are two non-equivalent optimal $`k`$-assignments $`\mu `$ and $`\nu `$. Assume further that

1. $`\mu `$ has a nonzero element in the last row and
2. $`\nu `$ does not use the last row.

Then there is a unique row $`s`$ such that

1. $`\nu `$ has a nonzero element in row $`s`$
2. There is a $`k`$-assignment $`\mu ^{}\mu `$ that does not use row $`s`$.

###### Proof.

Existence: Choose a $`k`$-assignment $`\mu ^{}\mu `$ which has as many matrix positions as possible in common with $`\nu `$. Since the last row is used by $`\mu ^{}`$ but not by $`\nu `$, there has to be a row $`s`$ which is used by $`\nu `$ but not by $`\mu ^{}`$. If $`\nu `$ has a zero element in row $`s`$, then $`\mu ^{}`$ must contain a zero element in the same column. By replacing this zero by the zero in row $`s`$, we would obtain a $`k`$-assignment $`\mu ^{\prime \prime }\mu ^{}`$ which has one more element in common with $`\nu `$, a contradiction. Therefore, the element in row $`s`$ which belongs to $`\nu `$ must be nonzero.

Uniqueness: Assume on the contrary that there are two different rows $`s`$ and $`t`$ that contain nonzero elements of $`\nu `$, say $`(s,c_1)`$ and $`(t,c_2)`$, and two equivalent $`k`$-assignments $`\mu \mu ^{}`$ such that $`\mu `$ does not use row $`s`$ and $`\mu ^{}`$ does not use row $`t`$. If $`\mu `$ would not use row $`t`$ either, then $`\mu \mathrm{}\nu `$ would have two distinct paths both with one end in a nonzero element of $`\nu `$, in the elements $`(s,c_1)`$ and $`(t,c_2)`$. This contradicts Lemma 3.2.

Hence there is a zero element $`(s,c_3)`$ of $`\mu ^{}`$, and similarly a zero element $`(t,c_4)`$ of $`\mu `$. Note that the symmetric difference $`\mu \mathrm{}\mu ^{}`$ contains only zeros and therefore it must consist of a number of $`(\mu ,\mu ^{})`$-alternating paths, all of even length. Let $`U`$ be the one that contains $`(t,c_4)`$. Unless $`U`$ ends at row $`s`$, $`\mu \mathrm{}U`$ will be a $`k`$-assignment equivalent to $`\mu `$ avoiding both rows $`s`$ and $`t`$. We have already seen that this is impossible.

We may therefore assume that $`\mu \mathrm{}\mu ^{}`$ consists of a single path from $`(t,c_4)`$ to $`(s,c_3)`$. Note that for every row $`r`$ which is used in this path, there is a $`k`$-assignment $`\mu ^{\prime \prime }\mu `$ that avoids row $`r`$. This is obtained by choosing the zeros from $`\mu `$ in the part of $`U`$ that goes towards $`(t,c_4)`$, and choosing them from $`\mu ^{}`$ in the part that goes towards $`(s,c_3)`$.

Let $`(m,d)`$ be the position in the last row used by $`\mu `$ and let $`L`$ be the $`(\mu ,\nu )`$-alternating path containing $`(m,d)`$. First note that if $`L`$ contains $`(s,c_1)`$, it has to end there and thus be of even length, which means that it has passed through $`(t,c_2)`$ first.

Case 1. $`L`$ intersects $`U`$. In this case, the first element of $`L`$ (starting from $`(m,d)`$) that belongs to $`U`$ must be an element of $`\mu `$. Suppose that this element is in row $`r`$. Then let $`\mu ^{\prime \prime }\mu `$ avoid row $`r`$. It follows that $`\mu ^{\prime \prime }\mathrm{}\nu `$ contains a path of even length starting at $`(m,d)`$, and another path containing the nonzero element $`(s,c_1)`$, contradicting Lemma 3.2.

Case 2. $`L`$ does not intersect $`U`$. Let $`S`$ be the path in $`\mu \mathrm{}\nu `$ containing $`(s,c_1)`$, $`(t,c_4)`$ and $`(t,c_2)`$ in this order.

We know that $`U`$ and $`S`$ intersect in $`(t,c_4)`$. Of the positions in $`US`$, let $`(r,c)`$ be the one which is closest to $`(s,c_3)`$ in $`U`$, and let $`\mu ^{\prime \prime }\mu `$ avoid row $`r`$. Then $`\mu ^{\prime \prime }\mathrm{}\nu `$ will consist of a cycle containing the nonzero position $`(s,c_1)`$, and two other paths, one containing $`(t,c_2)`$ and one containing $`(m,d)`$, contradicting Lemma 3.2. ∎

## 4. The probability that a row is used in the optimal assignment

Crucial for our proof is the following formula for the probability that a certain row without zeros is included in an optimal $`k`$-assignment. Let $`q(P)`$ denote this probability.

###### Theorem 4.1 (Row Inclusion Theorem).

Let $`P`$ be a standard RAP, and let $`r`$ be a row without zeros. The probability that some element in $`r`$ belongs to the optimal $`k`$-assignment is

(4)
$$q(P)=\frac{1}{m}\underset{i}{}\frac{\overline{d_{i,0}}(P)}{\left(\genfrac{}{}{0pt}{}{m1}{i}\right)},$$

where $`\overline{d_{i,0}}`$ is the number of partial $`k1`$-covers of $`i`$ rows not containing the row $`r`$.

To simplify notation we will assume that the nonzero row in the Row Inclusion Theorem is the *l*ast row in the matrix. A consequence of this formula is that the probability of using a certain row without zeros in an optimal $`k`$-assignment does not change if further zeros are introduced in a row which belongs to an optimal cover. This observation turns out to be sufficient for the proof of the formula.

###### Lemma 4.2.

Let $`P`$ be a standard RAP, and let $`r`$ be a row that belongs to an optimal cover. Let $`P^{}`$ be an RAP obtained from $`P`$ by inserting another zero in row $`r`$. Suppose there is a row without zeros in $`P`$. Then

$$q(P)=q(P^{}).$$

###### Proof.

Suppose that $`P`$ is a standard RAP where the first row belongs to an optimal cover, and that there is an element in the first row, say $`(1,1)`$, which is not zero. Suppose further that the last row contains no zeros. We want to show that if we replace the element $`(1,1)`$ by zero, the probability that the last row is used in the optimal $`k`$-assignment does not change.

Let $`\mathrm{\Omega }_P`$ be the probability space of all assignments of values to the random variables in $`P`$, that is, the space of all real nonnegative $`m`$ by $`n`$ matrices that have zeros in the positions where $`P`$ has zeros. If $`A\mathrm{\Omega }_P`$, we let $`q(A)`$ be 1 if the last row is used by an optimal $`k`$-assignment, and 0 otherwise. We let $`A_x`$, for nonnegative real $`x`$, denote the matrix obtained from $`A`$ by setting the entry in position $`(1,1)`$ to $`x`$.

We construct a measure preserving involution $`\phi `$ on $`\mathrm{\Omega }_P`$ with the property that (except possibly on a subset of probability zero) if $`A\mathrm{\Omega }_P`$ is a matrix where the last row changes between being used and not being used in the optimal $`k`$-assignment when $`A(1,1)`$ is set to zero, then in $`\phi (A)`$, the change goes the other way. In other words,

(5)
$$q(A)q(A_0)=q(\phi (A)_0)q(\phi (A)).$$

Let $`A\mathrm{\Omega }_P`$. If $`q(A)=q(A_0)`$, we let $`\phi (A)=A`$. Otherwise notice that if we slide the element in position $`(1,1)`$ continuously down to zero, there can be at most one point at which the location of the optimal $`k`$-assignment changes, and at this point, the element in position $`(1,1)`$ goes from not being used to being used.

At the point $`A_f`$ where the change occurs, the matrix is semi-generic with two non-equivalent optimal $`k`$-assignments. Let $`\mu `$ be the one that contains an element of the last row, and let $`\nu `$ be the one that doesn’t. By Lemma 3.4, there is a unique row $`s`$ such that $`\nu `$ contains a nonzero element in row $`s`$, and so that there is a $`k`$-assignment $`\mu ^{}\mu `$ that contains no element of row $`s`$. Notice that $`s`$ cannot be the first row, since the first row is used by every optimal $`k`$-assignment.

We let $`\phi (A)`$ be the matrix obtained from $`A`$ by swapping the elements in row $`s`$ with the corresponding elements (elements in the same column) in the last row, except in the columns where row $`s`$ has zeros.

In the analysis of this mapping, we introduce some auxiliary matrices. Let $`A_f^{}`$ be obtained from $`A_f`$ by setting the elements in the last row that are in columns where row $`s`$ has zeros, to zero. By Lemma 3.3, we can set these elements as small as we please without changing the location of the optimal assignments. Hence there will be no $`k`$-assignment of smaller cost than $`\mu `$ and $`\nu `$ in $`A_f^{}`$. Then since neither $`\mu `$ nor $`\nu `$ uses any zero element in row $`s`$, these can of course be increased without changing the optimality of $`\mu `$ and $`\nu `$. We let $`A_f^{\prime \prime }`$ be the matrix where the zero positions in row $`s`$ are changed to the values of the corresponding elements in the last row of $`A`$. If $`A`$ is generic, every optimal $`k`$-assignment in $`A_f^{\prime \prime }`$ is equivalent either to $`\mu `$ or to $`\nu `$.

Now $`\phi (A)_f`$ is obtained from $`A_f^{\prime \prime }`$ by swapping row $`s`$ with the last row. This means that in $`\phi (A)_f`$, there are exactly two equivalence classes of optimal $`k`$-assignments, one that includes position $`(1,1)`$ and one that doesn’t. But since the last row has been swapped with another row, we have, for $`xf`$,

$$q(\phi (A)_x)=1q(A_x),$$

which implies (5).

In $`\phi (A)_f`$, row $`s`$ has the property expressed in the conclusion of Lemma 3.4, but with the roles of $`\mu `$ and $`\nu `$ interchanged. It follows that $`\phi (\phi (A))=A`$, and in particular that $`\phi `$ is invertible. The mapping $`\phi `$ is piecewise linear, and on each piece, it is a permutation of variables. Hence it is measure preserving. This completes the proof. ∎

###### Proof of Theorem 4.1.

We now establish the Row Inclusion Theorem by an inductive argument. Let $`P`$ be a standard RAP, and suppose the last row contains no zero element. Notice that by König’s theorem (see \[LP86\]), the formula holds whenever $`P`$ has a $`k`$-assignment of only zeros. Suppose that the formula has been established for every standard RAP with fewer nonzero elements than $`P`$.

Case 1: Some optimal cover of $`P`$ contains a row $`r`$ with at least one nonzero element. Then every optimal $`k`$-assignment must use row $`r`$. Let $`P^{}`$ be as $`P`$ but with another zero in row $`r`$. By Theorem 4.2, the probability that the last row is used is the same in $`P`$ as in $`P^{}`$. Hence we only have to show that $`\overline{d_{i,0}}(P)=\overline{d_{i,0}}(P^{})`$ for $`i=0,\mathrm{},k1`$. We have to show that if a set of rows can be extended to a $`k1`$-cover of $`P`$, then it is possible to use row $`r`$ in this $`k1`$-cover. This will follow if we can show that row $`r`$ belongs to the optimal cover of the remaining zeros.

It suffices to show that if a row $`sr`$ is deleted from the matrix, row $`r`$ still belongs to the optimal cover of the remaining zeros. If the deletion of $`s`$ does not decrease the maximal number of independent zeros, this is obvious. Suppose therefore that the deletion of $`s`$ decreases the number of independent zeros. Then $`s`$ belongs to an optimal cover of $`P`$. Hence the row-maximal optimal cover of $`P`$ contains both $`r`$ and $`s`$. When row $`s`$ is deleted, the remaining rows and columns including $`r`$ will constitute an optimal cover of the remaining zeros.

Case 2: There is a row in $`P`$ with only zeros. Let $`P^{}`$ be obtained from $`P`$ by deleting this row and decreasing $`k`$ and $`m`$ by 1. The probability that the last row is used by an optimal assignment is clearly the same in $`P`$ as in $`P^{}`$. By induction, the probability that the last row is used in $`P^{}`$ is (with $`k=k(P),m=m(P)`$)

$$\frac{1}{m1}\underset{i=0}{\overset{k2}{}}\frac{\overline{d_{i,0}}(P^{})}{\left(\genfrac{}{}{0pt}{}{m2}{i}\right)}.$$

We have

$$\overline{d_{i,0}}(P)=\overline{d_{i,0}}(P^{})+\overline{d_{i1,0}}(P^{}),$$

for every $`i`$. Hence

(6)
$$\begin{array}{c}\frac{1}{m}\underset{i}{}\frac{\overline{d_{i,0}}(P)}{\left(\genfrac{}{}{0pt}{}{m1}{i}\right)}=\frac{1}{m}\underset{i=0}{\overset{k1}{}}\frac{\overline{d_{i,0}}(P^{})+\overline{d_{i1,0}}(P^{})}{\left(\genfrac{}{}{0pt}{}{m1}{i}\right)}\hfill \\ \hfill =\frac{1}{m}\underset{i}{}\frac{\overline{d_{i,0}}(P^{})}{\left(\genfrac{}{}{0pt}{}{m1}{i}\right)}+\frac{1}{m}\underset{i}{}\frac{\overline{d_{i,0}}(P^{})}{\left(\genfrac{}{}{0pt}{}{m1}{i+1}\right)}\\ \hfill =\underset{i}{}\overline{d_{i,0}}(P^{})\left(\frac{1}{m\left(\genfrac{}{}{0pt}{}{m1}{i}\right)}+\frac{1}{m\left(\genfrac{}{}{0pt}{}{m1}{i+1}\right)}\right)=\frac{1}{m1}\underset{i}{}\frac{\overline{d_{i.0}}(P^{})}{\left(\genfrac{}{}{0pt}{}{m2}{i}\right)}.\end{array}$$

Case 3: There is a unique optimal cover $`c`$ consisting of only columns. Then these columns will be used by the optimal $`k`$-assignment. Therefore the number of elements not covered by $`c`$ in the optimal $`k`$-assignment is independent of the random variables in the matrix. We condition on the position of the minimal element not covered by $`c`$. If in each case we subtract this minimum from all elements not covered by $`c`$, the same nonzero elements will be used in the optimal $`k`$-assignments. Here we are using Theorem 2.4 in the special case of no doubly covered elements. Let $`P^t`$ be the RAP obtained by conditioning on the minimal element in $`P`$ not covered by $`c`$ being in row $`t`$, and subtracting this minimum from all elements not covered by $`c`$. Then $`P^t`$ is a standard RAP, and the new zero occurring in row $`t`$ means that row $`t`$ belongs to an optimal cover of $`P^t`$. Hence in case the new zero is in the last row, that row must be used in the optimal $`k`$-assignment, while if it is not, we can find the probability that the last row is used by induction. We let $`P\backslash t`$ be the RAP obtained by deleting row $`t`$ from $`P`$ and decreasing $`k`$ by 1. Case 1 and 2 implies that the probability that the last row is used is the same in $`P^t`$ as in $`P\backslash t`$. It follows that the probability that the last row is used in the optimal $`k`$-assignment in $`P`$ is given by

(7)
$$\frac{1}{m}+\frac{1}{m(m1)}\underset{t}{}\underset{i}{}\frac{\overline{d_{i,0}}(P\backslash t)}{\left(\genfrac{}{}{0pt}{}{m2}{i}\right)}.$$

Here we have

$$\underset{t}{}\overline{d_{i,0}}(P\backslash t)=(i+1)\overline{d_{i+1,0}}(P),$$

since both sides are equal to the sum, taken over all $`t`$, of the number of partial $`k1`$-covers of $`P`$ that use row $`t`$ (and not the last row). It follows that (7) equals

(8)
$$\begin{array}{c}\frac{1}{m}+\frac{1}{m(m1)}\underset{i}{}\frac{(i+1)\overline{d_{i+1,0}}(P)}{\left(\genfrac{}{}{0pt}{}{m2}{i}\right)}\hfill \\ \hfill =\frac{1}{m}+\frac{1}{m}\underset{i}{}\frac{\overline{d_{i+1,0}}(P)}{\left(\genfrac{}{}{0pt}{}{m1}{i+1}\right)}=\frac{1}{m}\underset{i}{}\frac{\overline{d_{i,0}}(P)}{\left(\genfrac{}{}{0pt}{}{m1}{i}\right)}.\end{array}$$

This completes the proof of Theorem 4.1.

## 5. Proof of the cover formula for large $`m`$ or large $`n`$

The results of the last section will now enable us to prove that the cover formula (2) holds for standard RAP’s whenever either $`m`$ or $`n`$ is sufficiently large compared to $`k`$. We first prove that the cover formula is consistent with the Row Inclusion Theorem.

###### Theorem 5.1.

Let $`P`$ be a standard RAP where the last row contains no zeros, and let $`P_t`$ be obtained from $`P`$ by setting the element in column $`t`$ of the last row to zero. If the cover formula (2) holds for every $`P_t`$, then it holds for $`P`$.

###### Proof.

By Theorem 2.2 the probability that the element in column $`t`$ in the last row belongs to an optimal $`k`$-assignment in $`P`$ is equal to $`E(P)E(P_t)`$. Hence

$$nE(P)\underset{t}{}E(P_t)=\frac{1}{m}\underset{i}{}\frac{\overline{d_{i,0}}(P)}{\left(\genfrac{}{}{0pt}{}{m1}{i}\right)}.$$

Suppose that the cover formula holds for each $`P_t`$. Then

$$nE(P)=\frac{1}{mn}\underset{t}{}\underset{i,j}{}\frac{d_{i,j}(P_t)}{\left(\genfrac{}{}{0pt}{}{m1}{i}\right)\left(\genfrac{}{}{0pt}{}{n1}{j}\right)}+\frac{1}{m}\underset{i}{}\frac{\overline{d_{i,0}}(P)}{\left(\genfrac{}{}{0pt}{}{m1}{i}\right)}.$$

In order to prove that the cover formula holds for $`P`$, it is sufficient to prove that

(9)
$$\underset{i,j}{}\frac{d_{i,j}(P)}{\left(\genfrac{}{}{0pt}{}{m1}{i}\right)\left(\genfrac{}{}{0pt}{}{n1}{j}\right)}=\frac{1}{n}\underset{t}{}\underset{i,j}{}\frac{d_{i,j}(P_t)}{\left(\genfrac{}{}{0pt}{}{m1}{i}\right)\left(\genfrac{}{}{0pt}{}{n1}{j}\right)}+\underset{i}{}\frac{\overline{d_{i,0}}(P)}{\left(\genfrac{}{}{0pt}{}{m1}{i}\right)}.$$

If we write

$$\overline{d_{i,0}}(P)=\underset{j}{}\left(\frac{\overline{d_{i,j}}(P)}{\left(\genfrac{}{}{0pt}{}{n}{j}\right)}\frac{\overline{d_{i,j+1}}(P)}{\left(\genfrac{}{}{0pt}{}{n}{j+1}\right)}\right),$$

and fix $`i`$ and $`j`$, we see that (9) will follow from the identity

$$\frac{d_{i,j}(P)}{\left(\genfrac{}{}{0pt}{}{n1}{j}\right)}=\frac{\overline{d_{i,j}}(P)}{\left(\genfrac{}{}{0pt}{}{n}{j}\right)}\frac{\overline{d_{i,j+1}}(P)}{\left(\genfrac{}{}{0pt}{}{n}{j+1}\right)}+\frac{1}{n\left(\genfrac{}{}{0pt}{}{n1}{j}\right)}\underset{t}{}d_{i,j}(P_t).$$

Here, partial covers that contain the last row will contribute $`1/\left(\genfrac{}{}{0pt}{}{n1}{j}\right)`$ to both sides. If we let $`\overline{d_{i,j}}(P_t)`$ denote the number of partial $`k1`$-covers of $`P_t`$ with $`i`$ rows and $`j`$ columns that do not contain the last row, it only remains to show that

$$\frac{\overline{d_{i,j}}(P)}{\left(\genfrac{}{}{0pt}{}{n1}{j}\right)}=\frac{\overline{d_{i,j}}(P)}{\left(\genfrac{}{}{0pt}{}{n}{j}\right)}\frac{\overline{d_{i,j+1}}(P)}{\left(\genfrac{}{}{0pt}{}{n}{j+1}\right)}+\frac{1}{n\left(\genfrac{}{}{0pt}{}{n1}{j}\right)}\underset{t}{}\overline{d_{i,j}}(P_t),$$

or equivalently

$$j\overline{d_{i,j}}(P)+(j+1)\overline{d_{i,j+1}}(P)=\underset{t}{}\overline{d_{i,j}}(P_t).$$

Here the first term of the left hand side counts the partial $`k1`$-covers of $`P_t`$ that contain column $`t`$, while the second term counts those that don’t. Hence (9) holds. ∎

Next we show that the cover formula is consistent with removing a column which belongs to every $`k1`$-cover.

###### Theorem 5.2.

Let $`P`$ be a standard RAP, and suppose that the first column belongs to every $`k1`$-cover. Let $`P^{}`$ be the RAP obtained from $`P`$ by deleting the first column and decreasing $`k`$ and $`n`$ by 1. If the cover formula holds for $`P^{}`$, then it holds for $`P`$.

###### Proof.

By Theorem 2.5 $`E(P)=E(P^{})`$. We have (with $`n=n(P)`$ and $`k=k(P)`$)

(10)
$$\begin{array}{c}\frac{1}{mn}\underset{i,j}{}\frac{d_{i,j}(P)}{\left(\genfrac{}{}{0pt}{}{m1}{i}\right)\left(\genfrac{}{}{0pt}{}{n1}{j}\right)}=\frac{1}{mn}\underset{i,j}{}\frac{d_{i,j1}(P^{})+d_{i,j}(P^{})}{\left(\genfrac{}{}{0pt}{}{m1}{i}\right)\left(\genfrac{}{}{0pt}{}{n1}{j}\right)}\hfill \\ \hfill =\frac{1}{mn}\underset{i,j}{}\left(\frac{d_{i,j}(P^{})}{\left(\genfrac{}{}{0pt}{}{m1}{i}\right)\left(\genfrac{}{}{0pt}{}{n1}{j}\right)}+\frac{d_{i,j}(P^{})}{\left(\genfrac{}{}{0pt}{}{m1}{i}\right)\left(\genfrac{}{}{0pt}{}{n1}{j+1}\right)}\right)=\frac{1}{m(n1)}\underset{i,j}{}\frac{d_{i,j}(P^{})}{\left(\genfrac{}{}{0pt}{}{m1}{i}\right)\left(\genfrac{}{}{0pt}{}{n2}{j}\right)}\\ \hfill =E(P^{})=E(P).\end{array}$$

###### Theorem 5.3.

If $`P`$ is a standard RAP such that $`\mathrm{max}(m,n)>(k1)^2`$, then

$$E(P)=\frac{1}{mn}\underset{i,j}{}\frac{d_{i,j}(P)}{\left(\genfrac{}{}{0pt}{}{m1}{i}\right)\left(\genfrac{}{}{0pt}{}{n1}{j}\right)}.$$

###### Proof.

We show that if $`m>(k1)^2`$, then a standard RAP which contains a zero in each row must have a column which belongs to every $`k1`$-cover. The theorem then follows by induction on the number of nonzero elements of $`P`$. It suffices to consider the case that $`P`$ has exactly one zero in each row. Then the total number of zeros is greater than $`(k1)^2`$. Either there is a column with $`k`$ or more zeros, which must then belong to every $`k1`$-cover, or there is an independent set of $`k`$ zeros, and in this case there is no $`k1`$-cover.

By symmetry, the formula also holds if $`n>(k1)^2`$. ∎

## 6. Rationality of $`E(P)`$ as a function of $`n`$

We know from Section 5 that for fixed $`k`$, whenever $`m`$ or $`n`$ is large, the cover formula (2) for standard RAP’s holds. In order to prove that the formula holds for smaller values of $`m`$ and $`n`$, it is therefore sufficient to show that if $`k`$, $`m`$ and the zero positions are fixed, and we let $`P_n`$ be the standard RAP with $`n`$ columns, then there is a rational function in the variable $`n`$ that gives the value of $`E(P_n)`$ for every $`n`$ which is at least as large as $`k`$ and the number of columns with zeros. This rational function must then be equal to the one given by the cover formula. We prove this by induction over a class of RAP’s which includes not only standard RAP’s.

In an *exponential* RAP, all matrix elements are linear combinations with nonnegative rational coefficients of a set $`X_1,\mathrm{},X_p`$ of independent exponentially distributed random variables. If a variable $`X_i`$ in an exponential RAP has intensity 1, and occurs in one and only one matrix position, and this matrix entry is equal to $`X_i`$, then the variable is called a *standard variable*, and the matrix position where it occurs is called a *standard position*.

We introduce the concept of an *RAP-sequence*. The idea is to treat a set of similar RAP’s with different number of columns in a uniform way, in order to prove that there is a rational expression in the number of columns that gives the value of each RAP in the set.

###### Definition 6.1.

We say that a linear function $`f(x)=ax+b`$ in one variable is *$`k`$-positive*, iff $`f(x)>0`$ whenever $`xk`$, or equivalently, if $`a0`$ and $`b>ak`$.

Obviously, the sum of two or more $`k`$-positive functions is $`k`$-positive.

###### Definition 6.2.

An *RAP-sequence* is a sequence $`P_n`$ ($`nn_0`$) of exponential RAP’s satisfying the following:

1. Each $`P_n`$ is an exponential RAP in a set of variables $`X_{n,1},\mathrm{},X_{n,p(n)}`$.
2. The numbers $`m`$ and $`k`$ are uniform, that is $`m(P_i)=m(P_j)`$ and $`k(P_i)=k(P_j)`$ for all $`i`$ and $`j`$, while $`P_n`$ has $`n`$ columns.
3. In the first $`n_0`$ columns, the matrix elements of the various RAP’s differ only in that the first index of the variables is changed. In other words, the coefficient of $`X_{n,i}`$ in a matrix element in the first $`n_0`$ columns of $`P_n`$ is equal to the coefficient of $`X_{n_0,i}`$ in the same position in $`P_{n_0}`$.
4. Beyond column $`n_0`$, $`P_n`$ has only standard elements.

Moreover, we say that the RAP-sequence is *well-behaved* if

1. For each $`i`$ such that the variables $`X_{n,i}`$ occur in the first $`n_0`$ columns, there is a $`k`$-positive linear polynomial $`f_i(n)`$ such that the intensity of $`X_{n,i}`$ is $`f_i(n)`$.
2. Nonzero nonstandard elements occur only in columns which belong to the column-maximal optimal cover of the zeros, or equivalently, columns that intersect every maximal set of independent zeros.

###### Theorem 6.3 (Rationality Theorem).

Suppose $`P`$ is a well-behaved RAP-sequence. Then there is a rational function $`f(x)`$ in one variable such that

1. If $`x`$ is a zero of the denominator of $`f`$, then $`x<k`$.
2. $`E(P_n)=f(n)`$ for every $`nn_0(P)`$.

###### Definition 6.4.

We say that a linear combination $`u_1`$ of variables $`X_1,\mathrm{},X_p`$ is smaller than another if for each $`X_i`$, the coefficient is smaller in $`u_1`$ than in $`u_2`$. We say that $`u_1`$ and $`u_2`$ are *incomparable* if neither of them is smaller than the other, in other words, if each of them has a higher coefficient than the other for some variable.

We say that a linear combination is *potentially minimal* in a set of linear combinations, if it is smaller than or incomparable with every other.

Let $`c`$ be the row-maximal optimal cover of $`P`$. We prove the rationality theorem by induction on a number of parameters, in the following order:

1. The size of the largest independent set of zeros. An RAP-sequence is considered simpler if it has a larger set of independent zeros.
2. The number of rows in $`c`$. If the number of independent zeros are equal, the RAP-sequence with fewer rows belonging to the row-maximal optimal cover is simpler.
3. The set of potentially minimal nonstandard elements not covered by $`c`$. If 1. and 2. are equal, an RAP-sequence is considered simpler if fewer of the nonstandard elements not covered by $`c`$ are potentially minimal.
4. If 1–3 are equal, and there are two incomparable nonstandard elements not covered by $`c`$, then an RAP-sequence is simpler if there are fewer variables with different coefficients in the first two (in lexicographic order, say) incomparable potentially minimal nonstandard elements not covered by $`c`$.
5. If 1–3 are equal, and there is a minimal non-covered nonstandard element, then an RAP-sequence is simpler if the number of variables occurring in this element is smaller.

Let $`P`$ be a well-behaved RAP-sequence, and suppose that the rationality theorem holds for every simpler well-behaved RAP-sequence with the same values of $`m`$ and $`k`$. We may of course assume that $`|c|<k`$. Since $`P`$ is well-behaved, in each row all but at most $`k1`$ elements are standard and not covered by $`c`$.

We show that $`E(P)`$ can be expressed in terms of rational functions in $`n`$, and values of simpler well-behaved RAP-sequences.

Case 1: There are two or more non-covered incomparable nonstandard elements. Let $`u_1`$ and $`u_2`$ be the first two (in lexicographic order). We choose $`i`$ and $`j`$ such that the coefficient of $`X_{n,i}`$ is greater in $`u_1`$, and the coefficient of $`X_{n,j}`$ is greater in $`u_2`$. Let the coefficients be $`a_1`$, $`a_2`$, $`b_1`$ and $`b_2`$ so that $`u_1=a_1X_{n,i}+b_1X_{n,j}+\mathrm{}`$ and $`u_2=a_2X_{n,i}+b_2X_{n,j}+\mathrm{}`$.

Let $`Q`$ and $`R`$ be the RAP-sequences obtained by conditioning on $`(a_1a_2)X_{n,i}`$ being smaller or greater than $`(b_2b_1)X_{n,j}`$, respectively.

The intensities of $`(a_1a_2)X_{n,i}`$ and $`(b_2b_1)X_{n,j}`$ are $`f_i(n)/(a_1a_2)`$ and $`f_j(n)/(b_2b_1)`$ respectively. The probability of $`(a_1a_2)X_{n,i}`$ being smaller than $`(b_2b_1)X_{n,j}`$ is

$$\frac{f_i(n)/(a_1a_2)}{f_i(n)/(a_1a_2)+f_j(n)/(b_2b_1)},$$

and similarly, the probability of $`(b_2b_1)X_{n,j}`$ being smaller than $`(a_1a_2)X_{n,i}`$ is

$$\frac{f_j(n)/(b_2b_1)}{f_i(n)/(a_1a_2)+f_j(n)/(b_2b_1)}.$$

Therefore,

(11)
$$E(P_n)=\frac{f_i(n)E(Q_n)/(a_1a_2)+f_j(n)E(R_n)/(b_2b_1)}{f_i(n)/(a_1a_2)+f_j(n)/(b_2b_1)}.$$

We show that $`Q`$ and $`R`$ can be regarded as well-behaved RAP-sequences. If we condition on $`(a_1a_2)X_{n,i}`$ being smaller than $`(b_2b_1)X_{n,j}`$, then we can write

$$(a_1a_2)X_{n,i}=Y_n$$

and

$$(b_2b_1)X_{n,j}=Y_n+Z_n,$$

where $`Y_n`$ and $`Z_n`$ are independent exponentially distributed variables. The intensities are given by

$$I(Y_n)=\frac{f_i(n)}{a_1a_2}+\frac{f_j(n)}{b_2b_1}$$

and

$$I(Z_n)=\frac{f_j(n)}{b_2b_1},$$

both of which are $`k`$-positive.

When replacing $`X_{n,i}`$ and $`X_{n,j}`$ by the new variables $`Y_n`$ and $`Z_n`$, only the nonstandard elements are affected. Either $`u_1`$ is smaller than $`u_2`$ in $`Q`$, or at least the number of variables with distinct coefficients is smaller than in $`P`$, since $`u_1`$ and $`u_2`$ will get the same coefficient for $`Y_n`$. For any two matrix elements that satisfied $`u_3u_4`$ before we conditioned on which one of $`X_{n,i}`$ and $`X_{n,j}`$ is smallest, this inequality will still hold.

Hence $`Q`$, and similarly $`R`$, are simpler than $`P`$. By induction, it follows that (11) gives a rational expression for $`E(P)`$, where the denominator is nonzero for $`nk`$.

Case 2: There are no two non-covered incomparable nonstandard elements. Then either there is among the non-covered nonstandard elements a minimal one, or there are no non-covered nonstandard elements. We can treat these slightly different cases in the same way.

If there is a minimal non-covered nonstandard element, let $`aX_{n,i}`$ be a term occurring in this matrix entry. We let $`S_n`$ be a set of random variables consisting of $`aX_{n,i}`$ and all the non-covered standard variables (if there is no non-covered nonstandard element, we let $`S_n`$ consist only of the non-covered standard elements). There are at least $`nk+1`$ non-covered standard elements in each non-covered row. By grouping together the standard elements in each row, we can write the total intensity of $`S_n`$ in a uniform way as a sum of $`k`$-positive terms. Therefore the intensity of the minimum $`Y_n`$ of the terms in $`S_n`$ is $`k`$-positive as a function of $`n`$.

We condition on the minimal element in $`S_n`$. We can then replace the terms in $`S_n`$ by new variables $`Y_n`$ and $`Z_{n,i}`$ where $`Y_n=\mathrm{min}(S_n)`$ and $`Z_{n,i}`$ are the differences between the remaining terms in $`S_n`$ and the minimum.

Since all non-covered nonstandard elements contain the variable $`X_{n,i}`$ with a coefficient of at least $`a`$, we can by Theorem 2.4 subtract the minimum of $`S_n`$ from every non-covered element, and add it to the doubly covered elements.

We then get

(12)
$$E(P_n)=\frac{k|c|}{I(Y_n)}+\frac{1}{I(Y_n)}\underset{tS_n}{}I(t)E(Q_n(t)),$$

where $`Q_n(t)`$ are the new RAP’s obtained by conditioning on the term $`t`$ being smallest, and performing the change of variables and subtraction of the minimum, and $`I(t)/I(Y_n)`$ is the corresponding probability. In order to write this in a uniform way for the different values of $`n`$, we group together the cases where the minimum occurs in a particular row and beyond column $`n_0`$.

By permuting columns, we may assume that in those cases, the minimum always occurs in column $`n_0+1`$. The probability for each such case is $`(nn_0)/I(Y_n)`$. In this way, (12) will contain the same terms for each $`n>n_0`$.

If we let $`S_n^{}`$ be the subset of $`S_n`$ consisting of variables occurring in the first $`n_0`$ columns, and we let $`S_n^{\prime \prime }`$ be the set of standard variables in column $`n_0+1`$, we get

(13)
$$E(P_n)=\frac{k|c|}{I(Y_n)}+\frac{1}{I(Y_n)}\underset{tS_n^{}}{}I(t)E(Q_n(t))+\frac{nn_0}{I(Y_n)}\underset{tS_n^{\prime \prime }}{}I(t)E(Q_n(t)).$$

For $`n=n_0`$, the second sum will be empty. However, provided we can show that $`Q(t)`$ is a well-behaved RAP-sequence simpler than $`P`$, it will follow by induction that the denominator of $`E(Q(t))`$ does not vanish for $`nk`$. Hence the rational expression that occurs when multiplying the probability $`(nn_0)/I(S_n)`$ with the expression for $`E(Q(t))`$ will vanish for $`n=n_0`$. Therefore the rational expression that results will give the correct value of $`E(P)`$ also for $`n=n_0`$.

It remains to show that $`Q(t)`$ is a well-behaved RAP-sequence, and that it is simpler than $`P`$. We first show that $`Q(t)`$ is well-behaved. We have already seen that the intensities of the variables occurring in $`Q(t)`$ are $`k`$-positive. We therefore turn to the distribution of nonstandard elements. Possibly, there are some new nonstandard elements among the doubly covered elements. If there is no new zero among the non-covered elements, the columns of the doubly covered positions of course belong to the column-maximal optimal cover. Suppose that a new zero occurs among the non-covered elements. If the new zero is in a column that belongs to the column-maximal cover of $`P`$, then this is still the column-maximal cover of the new set of zeros. If the new zero is in a column that does not belong to the column-maximal cover, then since it is also in a row that does not belong to the row-maximal cover, there must be an independent set of zeros in $`Q(t)`$ which is larger than the largest independent set of zeros in $`P`$. Therefore, the column-maximal cover in $`P`$, extended with the column where the new zero has occurred, will be an optimal cover in $`Q(t)`$. Hence in any case, the columns of the column-maximal cover in $`P`$ belong to the column-maximal cover of $`Q(t)`$.

We now show that $`Q(t)`$ is simpler than $`P`$. If $`Q(t)`$ contains a new zero, then either it gives a larger independent set of zeros, or it has to be covered by a column in the row-maximal cover in $`Q(t)`$. In either case, $`Q(t)`$ is simpler than $`P`$. If no new zero occurs, this must be because the minimal term in $`S_n`$ was a term occurring in the minimal non-covered nonstandard element. In this case, the number of variables in this element will decrease, again making $`Q(t)`$ simpler than $`P`$.

Hence (12) gives a rational expression for $`E(P)`$ whose denominator is non-vanishing for $`nk`$.

This completes the proof of Theorem 6.3.

We are finally able to give a proof of Theorem 1.3, which we restate.

Theorem 1.3 *If $`P`$ is a standard RAP, then*

$$E(P)=\frac{1}{mn}\underset{i,j}{}\frac{d_{i,j}(P)}{\left(\genfrac{}{}{0pt}{}{m1}{i}\right)\left(\genfrac{}{}{0pt}{}{n1}{j}\right)}.$$

###### Proof.

A standard RAP with $`n_0`$ columns can be extended to a well-behaved RAP-sequence $`P_n`$ by inserting more columns without zeros. Hence Theorem 6.3 shows that there is a rational function $`f`$ giving the expected value of $`P_n`$. Since the cover coefficient $`d_{i,j}(P_n)`$ can be expressed uniformly as a polynomial in $`n`$, the cover formula gives a rational function in $`n`$ which takes the same values as $`f`$ on the infinitely many integers $`n>\mathrm{max}(n_0,(k1)^2)`$. Hence the cover formula must agree with $`f`$, and give the value of $`E(P)`$. ∎

## 7. Asymptotic results

Besides giving a new proof of Aldous’ $`\zeta (2)`$-limit theorem, Theorem 1.3 also makes the conjectured limits of \[LW00\] rigorous. It follows already from the conjecture of Coppersmith and Sorkin, now Theorem 1.2, that if $`\alpha ,\beta 1`$, then as $`k\mathrm{}`$, the value of the optimal $`k`$-assignment in a $`[\alpha k]`$ by $`[\beta k]`$-matrix of exp(1)-variables converges to

$$_\mathrm{\Delta }\frac{dxdy}{(\alpha x)(\beta y)},$$

where $`\mathrm{\Delta }`$ is the triangle with vertices in $`(0,0)`$, $`(1,0)`$ and $`(0,1)`$. For instance, when $`\alpha =1`$ and $`\beta =2`$ the limit is equal to

$$\frac{\pi ^2}{12}\frac{(\mathrm{log}2)^2}{2}.$$

It also follows that the value of a standard RAP with zeros in a region which is scaled up with $`k`$, $`m`$ and $`n`$ will converge to a similar integral. In particular it is shown in \[LW00\] that the cover formula implies that the limit value of a standard square RAP with zeros outside an inscribed circle is equal to $`\pi ^2/24`$.

A result presented as a conjecture in \[O92\], and proved in \[A01\], states that in the case $`k=m=n`$ with no zeros, as $`n\mathrm{}`$, the probability that the smallest element in a row belongs to the optimal assignment converges to $`1/2`$. We can now give an exact formula for this probability for finite $`k`$,

$$\frac{1}{2}+\frac{1}{2k}.$$

In the case $`k=m`$, the probability that the smallest element in a particular row belongs to the optimal assignment is equal to the probability that the smallest element in the entire matrix does. For arbitrary $`k`$, $`m`$, and $`n`$, the probability that the smallest element in the matrix belongs to the optimal $`k`$-assignment is

$$1\frac{k(k1)}{2mn}.$$

These theorems are obtained by combining the results of \[LW00\] and Theorem 1.3.


format: text

text:
# Optimal Techniques in Two-dimensional Spectroscopy: Background Subtraction for the 21st Century

## 1 Introduction

For more than 100 years, optical astronomers have employed long-slit spectrographs to study the internal physics of heavenly objects. Such data contain the target’s spectrum dispersed at every location along a single position angle on the sky (unless the target is an unresolved source). Unfortunately, at every location along the slit, one also collects photons from the night-sky emission. This background must be removed from the data in order to reveal the spectrum of the intended target. With the expansion of spectroscopy to wide fields-of-view, one can employ multi-slit aperture plates to collect spectra for many objects simultaneously, or observe with a very long slit. Over the past several years, important advances have been made in the art of background subtraction, though largely in the area of small-aperture spectroscopy, either with fibers or “micro-slits” (e.g., Kurtz & Mink, 2000; Glazebrook & Bland-Hawthorn, 2001; Viton & Milliard, 2003). However, or typical long- and multi-slit spectroscopy, the most common reduction procedures do not make optimal use of the data, and the resultant spectra with which one must perform one’s science suffer in quality compared to what is achievable with more modern techniques and modest computing power.

This paper outlines a new technique which makes optimal use of the data to accurately perform the subtraction of the unwanted background spectrum from two-dimensional spectra before any rebinning of the data is performed. In the procedure one makes full use of the spectrograph distortions to improve the sampling of the sky background spectrum. The quality of the background subtraction is completely insensitive to the magnitude of the distortions imposed by the spectrograph’s camera, to the severity of the curvature of the spectral lines that is caused by the dispersive element, or even to any tilting of individual slitlets in an aperture mask. Furthermore, by explicitly employing maps of the y-distortion and line curvature in a two-dimensional spectrum, the model two-dimensional background spectrum not only follows the same line curvature and y-distortion as the data, but the spectral features are sampled (pixelated) in exactly the same way as the features are sampled in the raw observations. As a result, when one subtracts the model from the data, there are absolutely no sharp residuals at the edges of night-sky emission lines. With more traditional methods, two-dimensional spectra must be rectified before one can perform the subtraction of the background; such rectification procedures introduce artifacts into the data, particularly when the sampling is poor. These artifacts can manifest themselves as sharp residuals at the edges of features with strong gradients, e.g., the night sky lines. Furthermore, in traditional methods, observers are required to identify (and remove) cosmic rays and other bad pixels before rebinning the data. With the method discussed in this paper, the sky subtraction is performed before the data have been rebinned, and, as a result, cosmic rays can be “cleaned” after the task of removing the sky background.

In the following sections, the method for fitting the two-dimensional night-sky background is described. Subsequently, this powerful new technique is applied to data riddled with cosmic rays and bright night-sky emission lines. The examples include data from three spectrographs, in which the data span a range of sampling from marginal (LRIS), to slightly under-sampled (NIRSPEC), to grossly under-sampled (MIKE). And finally the the advantages of this method over traditional methods will be summarized.

## 2 The Basic Procedure

Two-dimensional spectra, when ultimately imaged onto a detector, suffer from two problems that must be dealt with before or during the process of background subtraction: (1) the fact that the two-dimensional spectra are not aligned exactly along the rows (or columns) of the detector and are often curved with respect to the natural coordinate system of the detector (the y-distortion); and (2) the general tendency for dispersers to impose a wavelength-dependent line curvature onto the two-dimensional spectra (which may have already been tilted or curved if the slit has been so cut into the aperture plate). Furthermore, the coarse pixel sizes of modern detectors impose a third problem which normally limits the accuracy with which one could previously deal with the first two.

In order to discuss these issues and their resolution, we first define the image of two-dimensional spectroscopic data as $`P(x,y)`$, where $`(x,y)`$ are pixel coordinates in the system of the original image (e.g., a CCD frame). Because of distortions imposed by the optics, the spatial coordinate on the sky, $`y_t`$, for a given pixel $`(x,y)`$ is a non-linear function, $`y_t=Y(x,y)`$. Furthermore, the wavelength of light, $`\lambda `$, incident onto a pixel $`(x,y)`$ is also a non-linear function of image position. For the purposes of modeling the two-dimensional background spectrum, we are less concerned with the actual wavelength, $`\lambda `$, of incident light than we are with the fact that there exists a coordinate system, $`(x_r,y_t)`$, in which $`x_r`$ is a wavelength-dependent coordinate that is orthogonal to the spatial coordinate $`y_t`$. The transformation to this system is $`x_r=X(x,y_t)`$ such that the wavelength of light incident on a pixel can be written $`\lambda =L(x_r)`$, where $`x_r=X(x,Y(x,y))`$. Thus there exists a convenient coordinate system $`(x_r,y_t)`$ for which $`L/y_t=0`$. The transformations $`Y(x,y)`$ and $`X(x,y_t)`$ can be measured with great precision from comparison lamp spectra or from the night-sky features themselves (see, e.g., Kelson et al., 2000, for a description of a robust algorithm using FFTs and cross-correlations in order to make full use of the available data to map these distortions)111Note that that exact knowledge of the distortions does not free the observer from artifacts imposed on one’s data by the process of interpolation..

With knowledge of $`Y(x,y)`$ and $`X(x,y_t)`$, the traditional method of sky-subtraction required one to interpolate $`P(x,y)`$ onto a regular grid in $`(x_r,y_t)`$ before fitting the two-dimensional sky spectrum at every individual interval of $`\lambda `$ in the rebinned image. Unfortunately, the process of rebinning the data (1) introduces correlated noise; (2) smears cosmic rays and other bad pixels which might not have been flagged/cleaned beforehand; and (3) produces artifacts at the edges of sharp features. This last problem sometimes leads users to invoke high spatial orders in the fitting of the sky background in order to subtract such artifacts at the edges of sky lines. Furthermore, the rebinning of the data forces every sky spectrum used in the fit to have a common pixelation. This process limits one’s ability to accurately model the two-dimensional sky spectrum where gradients (e.g., $`P/x_r`$) can be large.

Instead of rebinning the data before performing the modeling of the two-dimensional sky spectrum, we propose that one should perform a least-squares fit to the sky spectrum using the original data-frame, in which the distortions and line curvature have not been removed. Using $`y_t=Y(x,y)`$ and $`x_r=X(x,y_t)`$ one should model the background spectrum in $`P(x,y)`$ as a function of the rectified coordinates $`(x_r,y_t)`$.

Each pixel in the original image represents an integral of the flux within a box the size of 1 pixel, but in the analysis each observed pixel’s location is assumed to be $`(x_r,y_t)`$, rather than $`(x,y)`$. In this way, each pixel samples the sky spectrum at a known sub-pixel position. Figure 1, in which a small section of an LRIS (Oke et al., 1995) spectrum is shown, demonstrates the utility of this change in coordinate systems. The left-hand panel shows a region around 5577Å in a short two-dimensional spectrum obtained using the 600 g/mm grating ($`1.28`$Å/pixel). The resolution is $`3.5`$ pixels (FWHM). Note in the image how pixelated the edges of the sky line are. Rebinning such data would lead to spatially periodic artifacts along the edge of the line. In the right-hand panel, the thick line shows the spectrum from one row, indicating how pixelated and poorly sampled the gradients in the line profile are in a given CCD row. The thin line shows the pixel values in the image $`P(x,y)`$ as a function of $`x_r=X(x,y_t)`$, revealing how well the CCD frame samples the sky spectrum. Such over-sampling is lost when one rebins the data.

In Figure 2, a larger section of the same LRIS frame is shown in the top panel. The middle panel shows $`P(x,y)`$ plotted as a function of $`X(x,y_t)`$. Note that many spikes in the data are visible. These discrepant points correspond to cosmic rays and other bad pixels. Fortunately, the sky spectrum is actually quite over-sampled in the rectified coordinate system, with many redundant measurements of the sky at nearly the same sub-pixel location.

Because of the redundancy in the sampling of the background spectrum, it becomes straightforward to identify the cosmic rays and bad pixels, even when they appear where gradients in the background spectrum are large. In the bottom panel, the thick line shows the spectrum sampled by a single row. The smooth line, shifted towards positive values, is a smoothed version of the data in the middle panel. The smoothing that was adopted was a running 30th-percentile within a window 30 data elements wide in $`x_r`$. In the bottom panel, we also show the running $`4\sigma `$ scatter (robustly determined) within the same sliding window. By comparing $`P(x,y)`$ with the smoothed spectrum in the middle panel, one can employ a simple $`\sigma `$-clipping algorithm to reject the cosmic rays and bad pixels, even where there are strong gradients in $`P(x,y)`$. Normally, the rejection/identification of cosmic rays or bad pixels involves the comparison of a pixel in an image with its nearest neighbors. However, because of the poor sampling of strong gradients in the background spectrum, one should compare a given pixel to those nearest only in $`x_r=X(x,Y(x,y))`$, i.e. sampled at the same sub-pixel interval. Thus, the smoothing is done on a copy of the array $`P(x,y)`$ that has been sorted in order of increasing $`x_r`$. Because the line curvature (plus the additional tilt of the slitlet) has improved the sampling of the night-sky spectrum, cosmic rays that fall where there are strong gradients in the sky spectrum can still be robustly identified, as other CCD rows have sampled the sky spectrum with similar sub-pixel sampling as at the location of a given cosmic ray.

While cosmic rays and other bad pixel values tend to be isolated sharp features, objects are generally extended in $`(x_r,y_t)`$. The top panel of Figure 3 shows a subsection of data for a slitlet containing a bright object. The middle panel shows every pixel $`P(x,y)`$ plotted, again, as a function of $`X(x,y_t)`$. Most of the pixels in the data only contain flux from the night-sky background, and these data points follow a locus in the figure with very small scatter. As in Figure 2, the cosmic rays are clearly visible above the well-defined background spectrum. Also visible, however, is a collection of points in the data which peak above the sky spectrum at regular intervals of $`X(x,y_t)`$. If we take the residuals of $`P(x,y)`$ from the robustly-smoothed version of the sky spectrum in the middle panel, and plot the statistical significance of those residuals against the spatial coordinate $`y_t=Y(x,y)`$, as in the bottom panel of the figure, the object pixels are clearly visible as a significant positive deviation from zero. By employing the $`\sigma `$-clipping described above for flagging cosmic rays, one also singles out those pixels that are significantly contaminated by flux from the object. In general a choice of clipping at 3 or $`4\sigma `$ is very effective at rejecting cosmic rays from the fit to the sky background, and it also rejects objects which could seriously affect the fit to the sky. Faint objects, which do not peak above the adopted threshold tend not to affect the modeling of the sky and most such objects tend to cover too small a spatial area to adversely affect the fit anyway. While a clipping method works sufficiently well for most applications, one can straightforwardly implement an input set of sky apertures, outside of which pixels are simply ignored. Sky apertures may be specifically useful in very short slitlets, where the the sky covers $`\stackrel{<}{}10\%`$ of the spatial extent of the slit.

Once the discrepant pixels are rejected, a bivariate B-spline (e.g., de Boor, 1978; Dierckx, 1993) is constructed as an approximation to the remaining remaining data points as a function of their positions. However, the pixel values are not considered to represent a bivariate function of $`(x,y)`$ (the original CCD coordinates) but the data are treated as a bivariate function of $`(x_r,y_t)=(X(x,Y(x,y)),Y(x,y))`$. In the author’s implementation of the method, the DIERCKX surface- and curve-fitting library available from NETLIB was used. This library allows one to weight each pixel during the fit for the B-spline coefficients, and these weights were set equal to the inverse of the expected noise.

When fitting for the B-spline representation of a bivariate dataset, the smoothness of the model is set by the density of knots in the two cardinal directions. In the simplest construction of the two-dimensional sky background, the knot locations, $`t_x`$ and $`t_y`$, are chosen with a high density in $`x_r`$ such that $`t_x=\{\mathrm{min}(x_r),\mathrm{min}(x_r)+1,\mathrm{},\mathrm{max}(x_r)\}`$ with intervals of $`\mathrm{\Delta }x_r=1`$. The knots in $`y_t`$ are chosen to be $`t_y=\{\mathrm{min}(y_t),\mathrm{max}(y_t)\}`$, where the minima and maxima in $`x_r`$ and $`y_t`$ are derived for the slit being analyzed. In this way, the B-spline is non-parametric function in $`x_r`$ and a polynomial of order $`k_y`$ in $`y_t`$. The choice of $`k_y\{1,3,5,\mathrm{}\}`$ sets the order of the spatial variation in the sky spectrum at a fixed $`x_r`$ (e.g., fixed $`\lambda `$).

With the default choice of knots in $`x_r`$ described in the previous paragraph, the fit to the data very nearly approximates an interpolating spline along the wavelength-dependent coordinate. However, the B-spline was generated using data with finer sampling than that available in a single CCD row. As a result, the spline is a smooth representation of the sky spectrum at those locations $`(x_r,y_t)`$ in the original data frame. One can reduce the number of knots in $`t_x`$ to impose greater smoothness on the model sky spectrum in those ranges of $`x_r`$ with little structure, while leaving a higher density of knots near bright night-sky emission lines in order to better match the gradients there. While such optimization of the knots in $`t_x`$ can improve the background modeling with particularly noisy data, the selection of $`t_x`$ described above is generally sufficient. One can also insert more knots in $`t_y`$ to model more complicated spatial variation in the sky (at fixed $`\lambda `$). With clever placement of knots in the spatial direction, the model can better map high-frequency spatial variations such as in data containing residual fringes in the very red portions of a spectrum. If the spatial variation in the sky spectrum is expected to be negligible, such as in very short slits, a one-dimensional B-spline can be fit to the values of $`P(x,y)`$ as a univariate function of $`x_r`$.222One optimization for $`t_x`$ in the bivariate method includes first finding the optimal univariate B-spline representation of the sky spectrum, with the resulting optimal knot locations subsequently used in the bivariate fit to the data.

Regardless of the complexity of the knot placement, the B-spline representation of the data is well-behaved because every pixel location in the original image, $`(x,y)`$ has a correspondingly unique, rectified coordinate $`(x_r,y_t)`$. As a result the bivariate fit is well-defined and computationally inexpensive to construct. Once the B-spline representation has been computed, it can be quickly evaluated at every $`(x_r,y_t)`$ location in the spectrum being analyzed. In this way, the two-dimensional background spectrum is only evaluated at the locations where the original pixels exist. Therefore no interpolation at sub-pixel locations (in $`x,y`$) occurs in the original data frame. As a result, there are no artifacts at strong gradients.

Some examples of the application of this method of optimal background subtraction are shown in the next section. The examples shown below used the simplest knot selection in the above discussion. No optimization of knot location was implemented. Finally, the examples were generated using $`k_x=3`$ and $`k_y=1`$.

## 3 Examples

### 3.1 LRIS

The top panel of Figure 4 shows a section of two-dimensional spectrum from one strongly tilted slitlet. The wavelength range covers from blueward of Na I to approximately 6300Å. The second panel shows the model sky spectrum, derived using the new method. The third panel shows the difference between the two (i.e., the background-subtracted spectrum). Note how well the night-sky emission lines are subtracted, with no residual systematic structure. The bottom panel shows an $`rms`$-smoothed version of the background-subtracted image, normalized by the expected noise. This representation of the data illustrates that there is no additional noise at the locations of the sharp sky lines.

Figure 5 follows a similar format as Figure 4 but for a different slitlet, and for a wavelength range 7000Å $`\stackrel{<}{}\lambda \stackrel{<}{}`$ 7700Å. As in the previous figure, the background subtraction is free of any residual systematic structure, and the noise in the background-subtracted image is as expected.

Figure 6 shows a section of data with a wavelength range 7200Å $`\stackrel{<}{}\lambda \stackrel{<}{}`$ 7600Å. This section has a very strong cosmic ray at $`(x,y)(1220,15)`$ that sits along a large portion of a sky line. Note how clean the sky-subtracted image is, again, with no systematic residuals at the locations of bright lines.

Figure 7 shows a section of data with a wavelength range 5500Å $`\stackrel{<}{}\lambda \stackrel{<}{}`$ 6500Å. These data were taken from a slitlet at the edge of the LRIS field, at which the $`y`$-distortion, $`Y(x,y)`$, has a very strong derivative with respect to $`x`$. Note how cleanly the bright sky lines 5577Å, 5890Å, 5896Å, 6300Å, etc., are modeled and subtracted. As in the previous examples, there are no systematic residuals at the locations of bright night-sky features.

### 3.2 NIRSPEC

Figure 8 shows a two-dimensional $`H`$-band spectrum obtained with NIRSPEC (McLean et al., 1998). The dispersion is approximately 2.8Å/pixel. The distortions in the data are large, and the sampling is poor. Small sections of these data are also shown in Figures 9 and 10. In this example, the knots in the dispersion direction were located at half-pixel intervals. This choice was motivated by the strong distortions and length of the slit. Together these allow for higher resolution in the reconstructed background spectrum.

Even in coarsely sampled near-IR spectra, the method provides a two-dimensional model of the background with identical sampling as in the original data. Thus, the sky background subtracts cleanly, with no systematic residuals at the locations of the night sky emission lines. When the sampling is this coarse, the observer should avoid interpolating sharp features, and observers now can choose whether to rebin the sharp features in the raw data or to rebin the sky-subtracted frame. Figures 11 and 12 employ these NIRSPEC data to directly compare this new methodology for sky-subtraction with traditional procedures. In the figures, we show the NIRSPEC data in two states: (1) sky-subtracted first and rectified second; and (2) rectified first and sky-subtracted second. In the latter form, strong, periodic artifacts are plainly visible at the locations of the night sky emission lines, where the under-sampled edges were rebinned. In the former, where the task of sky-subtraction was performed first, no sharp gradients remain in the data to introduce artifacts into the rectified frame. While this example shows a comparison of rectified two-dimensional sky-subtracted data, some readers may choose not to rebin their two-dimensional spectra at all, and opt for extracting one-dimensional spectra directly from the unrectified background-subtracted frames (using knowledge of the two-dimensional wavelength solution).

### 3.3 MIKE

Figures 1315 show data from a one hour integration with the red side of the MIKE echelle spectrograph on the Clay 6.5m telescope at Magellan. Because the optical path in MIKE (Bernstein et al., 2003) involves using one prism as the primary disperser and the cross-disperser, the orders have large curvature and the spectra themselves show large line curvature. This exposure covers order #62 (bottom, central wavelength $`5540`$Å) through order #33 (top, central wavelength $`10300`$Å). The data had been binned $`2\times 2`$ while reading the CCD to avoid being read-noise dominated, even with the one hour integration. The binning increases the fraction of the image contaminated by cosmic rays by a factor of four, and reduces the spatial extent of the slit to about 20 pixels. Many users also use $`1\times 3`$ and even $`2\times 4`$ binning with MIKE, and as a result, most data obtained with MIKE will be severely under-sampled, and the binning leads to higher contamination rates with cosmic rays.

In Figure 15 we show a direct comparison of a portion of the rectified sky-subtracted MIKE data (similar to the comparison for NIRSPEC data in Figure 11). While at such high resolution the night sky lines do not affect much of the data, one may still wish to recover accurate estimates of the object flux at the locations of the sky lines. In the right-hand panel of Figure 15 there are clear artifacts at the locations of the sky-lines and such artifacts render the science data useless at those wavelengths. In the left-hand panel, the object spectrum is not contaminated by the artifacts introduced by the rebinning of night sky lines.

The examples shown in this section illustrate a range of results from the modeling procedure outlined in this paper. In all cases, when the distortions and line curvature are accurately known, the sky background shows very little spatial variation when fit in the rectified coordinate system (see §2). The need for a high spatial order is only required when using traditional methods, and artifacts at the edges of sky lines are not well modeled with low-order polynomials. Because of this, observers can also apply this method to short slitlets using one-sided sky, restricting the B-spline to one dimension (wavelength) or perhaps only limiting the spatial order to $`k_y=1`$. In such cases, the only factor limiting the accuracy of the sky subtraction would be flat-fielding and other techniques may prove more valuable (e.g., Glazebrook & Bland-Hawthorn, 2001).

## 4 Summary

This paper describes a new method for modeling the background in two-dimensional slit spectroscopy.333An implementation of the method is available as part of a suite of Python reduction scripts written by the first author at http://www.ociw.edu/$``$kelson. Though no documentation currently exists, readers may download at their discretion. This method makes full use of the data in its “raw” state — before the data have been rectified (rebinned). The distortions inherent in the data provide the means by which the background spectrum can be accurately modeled at repeated sub-pixel locations. Because of the improved sampling of the background spectrum in the original data compared to a rectified image, one can also robustly reject pixels contaminated with cosmic rays and bright objects. Because of these advantages over traditional background-subtraction techniques, the task of subtracting the two-dimensional background spectrum should now be performed as one of the first steps in any pipeline of spectroscopic reduction, even before cosmic rays and other bad pixels have been “cleaned”.

The subtraction of the background can now be performed on raw spectra, riddled with cosmic rays. As a result, the procedure can be adopted as one of the first steps in a pipeline for spectral reduction, and not as one of the last. The task is computationally inexpensive and can allow observers to quickly analyze incoming data at the telescope, assuming the distortion maps have been characterized from calibration data during the previous afternoon, or are known a priori from an optical model of the instrument. After the background has been subtracted the cosmic rays can be flagged trivially, without the need for complicated cosmic ray identification procedures.

Because the algorithm is utilizing the full sampling of the background spectrum produced by the distortions in the raw data, the method is making optimal use of the available data. While this procedure has been presented in the context of taking full advantage of the optical distortions and spectral line curvature to recover better sampling of areas with strong gradients in the data, the method works just as well on data for which the line curvature is negligible. This recovery of the two-dimensional background spectrum leaves no systematic residuals at the locations of strong gradients simply because the model is never interpolated at locations where the data do not exist. The model is always sampled where the original data exist. As a result, the sharp gradients, such as those at the edges of bright sky emission lines, are never rebinned and ringing does not occur.

For spectrograph/detector combinations that over-sample the data (e.g., DEIMOS and the IMACS long camera), traditional methods may be deemed satisfactory by many observers. However, in near-IR or high-resolution echelle spectroscopy, the data may be binned and/or under-sampled in order to reduce the importance of read noise. Regardless of the source of one’s spectroscopic data, observers now can choose whether to rebin their data with or without the sharp night-sky features present. Observers may now choose the latter, as the rebinning of coarsely sampled features produces unwanted artifacts that prevent the accurate modeling of the background, and prevent the accurate recovery of an object’s flux at every observed wavelength.

Nevertheless, for observers who wish to proceed from the telescope to extracted spectra in the fewest number of steps possible, without the additional correlated noise introduced by rebinning ones data, and without any artifacts caused by rebinning any strong gradients, computational machinery has now caught up to your demands.

The author would like to thank the staff at the Carnegie Observatories for their support, and in particular P. Martini and M. Rauch for the use of their NIRSPEC and MIKE data. Furthermore, the anonymous referee is acknowledged for ensuring that this work is received by a broader audience.