Inverting the matrix $A[r,c] = \frac{1}{r + c}$

 Suppose $A$ is an $n \times n$  matrix defined by

$$A[r,c] = \frac{1}{r+c}, \quad 1 \le r \le n, 1 \le c \le n$$

i.e the entry in the $r^{th}$ row and $c^{th}$ column is $\frac{1}{r+c}$.

i.e

$$ A = \begin{bmatrix} \frac{1}{2} & \frac{1}{3} & \dots & \frac{1}{n+1}\\ \vdots & \vdots & \dots & \vdots\\ \frac{1}{n+1} & \frac{1}{n+2} & \dots & \frac{1}{2n} \end{bmatrix}$$

Determine the inverse of $A$.

Scroll down for a (somewhat suprising) solution.

.

.

.

.

.

Solution:

[In the below $m$ and $i$ will denote a row index and $k$ and $j$ a column index (instead of $r$ and $c$).]

Fix a $k$ ($ 1 \le k \le n$).

The $k^{th}$ column of the inverse is the solution $(a_1, a_2, \dots, a_n)$ to the system of $n$ equations

$$ \sum_{j = 1}^{n} \frac{a_j}{m + j} = \delta_{mk} ,\quad 1 \le m \le n$$

where $$\delta_{mk} = \begin{cases} 0 & m \neq k \\ 1 & m = k \end{cases}$$

$\delta_{mk}$ are the entries of the identity matrix, basically and the LHS in the above is the dot product of the $m^{th}$ row of $A$ and $k^{th}$ column of $A^{-1}$.

Now consider 

$$f(x) = \sum_{j=1}^{n} \frac{a_j}{x + j} = \frac{P(x)}{Q(x)}$$

where $P(x)$ is an $(n-1)^{th}$ degree polynomial and $Q(x) = (x+1)(x+2)\dots (x+n)$.

The $n$ equations (with $\delta_{mk}$) basically mean that $1, 2, \dots, k-1, k+1, \dots, n$ are roots of $P(x)$.

Let $R(x) = (x-1)(x-2)\dots(x-n)$. Then we must have that

$$P(x) = A\frac{R(x)}{x - k}$$ 

for some constant $A$.

We also have $P(k) = Q(k)$  (because $f(k)  = \delta_{kk} = 1$) and thus $$Q(k) = P(k) =  \lim_{x \to k} A\frac{R(x)}{x-k} = A R'(k)$$

Which implies $$A = \frac{Q(k)}{R'(k)}$$

$R'(x)$ being the derivative of $R(x)$.

Thus we get the formula

$$ f(x) = \sum_{j=1}^{n} \frac{a_j}{x+j} = \frac{Q(k) R(x)}{R'(k)(x-k)Q(x)}$$ 

Now just looking at

$$ f(x) = \sum_{j=1}^{n} \frac{a_j}{x+j}$$

we get

 $$a_m = \lim_{x \to -m} f(x) (x+m)$$

i.e

$$a_m = \frac{Q(k)}{R'(k)} \lim_{x \to -m} \frac{R(x) (x+m)}{(x-k) Q(x)}$$

Now $$\lim_{x \to -m} \frac{x+m}{Q(x)} = \frac{1}{Q'(-m)}$$

Thus we get

$$a_m = -\frac{Q(k) R(-m)}{(m+k) R'(k) Q'(-m)}$$

Thus the inverse matrix $A^{-1}$ satisfies

$$A^{-1}[m, k] = -\frac{Q(k) R(-m)}{(m+k) R'(k) Q'(-m)}$$

Remember, $Q(x) = (x+1)(x+2)\dots (x+n)$ and $R(x) = (x-1)(x-2)\dots(x-n)$.

It turns out that we can rewrite the above as products of binomial coefficients, and they are all integers!

$\textbf{Additional Remarks:}$ Apparently our $A$ is a special case of a Cauchy Matrix, where $A[i,j] = \frac{1}{x_i - y_j}$ and a method similar to above works for the general case too and was solved in 1959!