Determinant of a generalized Pascal matrix Let $M$ denote the infinite matrix defined recursively by 
$$
M_{ij} = 
\begin{cases}
  1, &  \text{if } i=1 \text{ and } j=1; \\
  aM_{i-1,j}+bM_{i,j-1}+cM_{i-1,j-1}, & \mbox{otherwise}.\\ \end{cases}
$$
($M_{i,0}$ and $M_{0,j}$ are both defined to be $0$.) 
(Added: I just discovered that the numbers in the $M$ matrix are called weighted Delannoy numbers.)
Let $M_n$ denote the $n \times n$ upper-left submatrix of $M$. 
For example, with $a = b = c = 1$, 
$M_1 = \begin{bmatrix} 1 \end{bmatrix}$, $M_2= \begin{bmatrix} 1 & 1 \\ 1 & 3 \end{bmatrix}$, and $M_3 = \begin{bmatrix} 1 & 1 & 1 \\ 1 & 3 & 5 \\ 1 & 5 & 13 \end{bmatrix}$, and $M_4 = \begin{bmatrix}
    1 & 1 & 1 & 1 \\
    1 & 3 & 5 & 7 \\
    1 & 5 & 13 & 25 \\
    1 & 7 & 25 & 63
   \end{bmatrix}.$  
A few years ago one of my students proved, by induction, that $$\det M_n = (ab+c)^{n(n-1)/2}.$$
My question is

Is there a noninductive proof that $\det M_n = (ab+c)^{n(n-1)/2}$ that gives more insight into why the determinant works out so nicely?

For example, when $a = b = 1$, $c = 0$, $M$ is the symmetric Pascal matrix.  I've seen more than one way to prove that $\det M_n = 1$ in this case.  For example, Edelman and Strang give four proofs of an LU-decomposition that does it.  I also once saw, at a conference, a combinatorial proof using the interpretation of the determinant in terms of nonintersecting paths in a directed graph.  (I think the talk was given by Art Benjamin, but it was several years ago, and I may be misremembering.)  So I know that there are some nice proofs in the special case of the Pascal matrix.  But what about the general case?
 A: Inspired by the solution of Sivaram Ambikasaran, I would just like to furnish some concision. Given
$$
M_{i,j} = 
\begin{cases}
  a^{i-1}b^{j-1} &  \text{if } i=1 \text{ or } j=1, \\
  aM_{i-1,j}+bM_{i,j-1}+cM_{i-1,j-1} & \mbox{if }i>1 \text{ and } j>1,\\ \end{cases}
$$
one may (inspired by what is said in the Edelman-Strang paper) start subtracting from each row (except the first) $a$ times the previous row. This is to be done simultaneously, or working from the bottom up, so that it is the original value of the row that is subtracted. This amounts to left-multiplying by the elimination matrix $E_{-a}$ with diagonal entries $1$ and subdiagonal entries $-a$. One obtains a matrix $M'$ given by 
$$
M'_{i,j} = 
\begin{cases}
  b^{j-1} & \text{if } i=1, \\
  M_{i,j}-aM_{i-1,j} & \mbox{if }i>1.\\ \end{cases}
$$
The latter expression is $0$ if $j=1$, and otherwise
$$
  M'_{i,j} = bM_{i,j-1}+cM_{i-1,j-1}\qquad\text{if }i,j>1
$$
using the recursion defining $M$. Now proceed similarly by columns, right-multiplying by the transpose $E_{-b}^\top$ of $E_{-b}$, giving a new matrix $M''$ given by
$$
M''_{i,j} = 
\begin{cases}
  \delta_{i,j} & \text{if } j=1, \\
  M'_{i,j}-bM'_{i,j-1} & \mbox{if }j>1.\\ \end{cases}
$$
The latter expression is $0$ if $i=1$, and otherwise
$$\begin{align}
  M''_{i,j} &= bM_{i,j-1}+cM_{i-1,j-1}-b(M_{i,j-1}-aM_{i-1,j-1})\\
           & =(ab+c)M_{i-1,j-1}
\end{align}\qquad\text{if }i,j>1,
$$
using the two expressions for $M'_{i,j}$ given above. In other words one has in block form
$$
  M''=\begin{pmatrix}1&0\\ 0&(ab+c)M_{(n-1)}\end{pmatrix}
$$
where $M_{(n-1)}$ is the $(n-1)\times(n-1)$ top-left submatrix of $M$, which is its counterpart of size one less. Therefore
$$
  \det M_{(n)} = \det M = \det M''= (ab+c)^{n-1} \det M_{(n-1)}
$$
from which it follows by induction that $\det M=(ab+c)^{\tbinom{n}2}$.
The scalar multplication of $M_{(n-1)}$ by $ab+c$ can of course be realised as the multiplication by a multiple of the identity matrix, which commutes with with the matrices obtained in a recursive decomposition of $M_{(n-1)}$. It is then easy to conclude that one has a decomposition
$$
  M_{(n)} = L_{(n)} D_{(n)} U_{(n)}
$$
where $D_{(n)}$ is diagonal with entries $(D_{(n)})_{i,i}=(ab+c)^{i-1}$ and $L_{(n)}$ and $U_{(n)}$ are unitriangular and therefore have determinant $1$. Their explicit form can be obtained from the recurrences
$$
  L_{(n)} = (E_{-a})^{-1}\cdot \begin{pmatrix}1&0\\ 0&L_{(n-1)}\end{pmatrix},
  \qquad
  U_{(n)} = \begin{pmatrix}1&0\\ 0&U_{(n-1)}\end{pmatrix}\cdot(E_{-b}^\top)^{-1},
$$
which can easily found to be solved by $(L_{(n)})_{i,j}=\binom{i-1}{j-1}a^{i-j}$ and $(U_{(n)})_{i,j}=\binom{j-1}{i-1}b^{j-i}$, either using $((E_{-a})^{-1})_{i,j}=a^{i-j}$ or by solving first $((L_{(n)})^{-1})_{i,j}=\binom{i-1}{j-1}(-a)^{i-j}$. 
A: This is not an answer but rather a comment of Sivaram's answer.
[update] I've to correct the numerical examples because I used a slightly distorted version of M. But this does not affect the basic idea of using the LDU-decompostion instead of the LU-decomposition. Finally, this is now not much more than a slight optimization of Sivaram's ansatz
If I decompose into $\small L \cdot D \cdot U $ instead, where D is diagonal, then the pattern of the construction of the factors L and U is more "simple"(contains only a resp. b) but their diagonals are 1 and so can be neglected for the determinant. So we need only the determinant of D. Then we find, that D can be described by
$\qquad \small D=[1,c+ab, (c+ab)^2, (c+ab)^3,\ldots ,(c+ab)^{n-1} ] = dV(c+ab,n) $ 
and the determinant for matrix-size n x n is then the product of these terms up to n-1 by an obvious simple expression. Again the found pattern in the LDU-decomposition must be proven, but the approach might possibly allow a shorter path...   
[added] : The actual LDU-decomposition can be written as                    
$\qquad \small M = P^{\quad a} \cdot \quad dV(c+ab) \cdot (P^{\quad b} \sim )$     
where P is the lower triangular Pascalmatrix,  "~" means the transpose and dV(x) is the diagonalmatrix containing the consecutive powers of its argument x beginning at $\small x^0 $        
(remark: I also prefer this LDU often because it avoids introducing sqaures and/or squareroots, like for instance does the cholesky-decomposition for the symmetric case)

in case your software has a LDU-decomposition not available, here is some code usable for Pari/GP    

 LDU(Y) = local(dim=#Y, D, MR, ML, dx); \\ Y must be square, 
       \\ no errorchecks in     demo-documentation
 D=matrix(dim,dim); MR=matid(dim); ML=matid(dim);
 for(p=1,dim, 
   D[p,p]=(dx=Y[p,p]);
   for(c=p+1,dim,
        MR[p,c]=if(dx==0,0,Y[p,c]/dx)
       );
   for(r=p+1,dim, 
        ML[r,p]=if(dx==0,0,Y[r,p]/dx)
      );
   for(r=p+1,dim,
      for(c=p+1,dim,
            Y[r,c]-=ML[r,p]*dx*MR[p,c]
      ));
  );
return([ML,D,MR]);

The matrix L has the form [now corrected!] $\small L = P^a $   
$\small \begin{array} {lllll}
 1 & . & . & . & . & . \\
 a & 1 & . & . & . & . \\
 a^2 & 2a & 1 & . & . & . \\
 a^3 & 3a^2 & 3a & 1 & . & . \\
 a^4 & 4a^3 & 6a^2 & 4a & 1 & . \\
 a^5 & 5a^4 & 10a^3 & 10a^2 & 5a & 1
 \end{array} $
while U is simply the transposed and a replaced by b.     
 
[update]
Here is the Pari/GP code to generate the matrix M (corrected) and the call for the LDU-decomposition:

M = matrix(6,6,r,c); for(r=1,6,M[r,1]='a^(r-1)); for(c=1,6,M[1,c]='b^(r-1));
    for(r=2,rows(M),for(c=2,cols(M),M[r,c]='a*M[r-1,c]+'b*M[r,c-1]+'c*M[r-1,c-1]))
M_ldu = LDU(M); L=M_ldu[1];D=diag(M_ldu[2]);U=M_ldu[3];
print (D)  \\ check the diagonal-component
%1164 = [1,
         c + a*b ,
         c^2 + 2*b*a*c + b^2*a^2,
         ... (snipped) ] ~

A: The identity follows almost immediately from the determinant formula for counting non-intersecting paths:


*

*on one hand, the determinant counts the number of non-intersecting paths from points $a_i=(-i,0)$ to points $b_i=(0,i)$ on the lattice with $a$ edges going right, $b$ up, and $c$ up-right (for each vertex);

*on the other hand, to choose non-intersecting paths means to choose path from each of $n(n-1)/2$ "internal" vertices to it's upper-right neighbor — which can be done in exactly $ab+c$ ways: one can go either up-then-right or diagonally (for example, when $c=0$ there is, obviously, just one "geometric" path with some weight). Upd: that's not quite true; I'll try to fix it.
