  This article deals with the study of the following singular quasilinear
equation:
  \begin{equation*}
  (P) \left\{ \ -\Delta_{p}u -\Delta_{q}u = f(x) u^{-\delta},\; u>0 \text{ in
}\; \Om; \; u=0 \text{ on } \pa\Om,
  \right.
  \end{equation*}
  where $\Om$ is a bounded domain in $\mathbb{R}^N$ with $C^2$ boundary
$\pa\Om$, $1< q< p<\infty$, $\de>0$ and $f\in L^\infty_{loc}(\Om)$ is a
non-negative function which behaves like $\textnormal{dist}(x,\pa\Om)^{-\ba},$
$\ba\ge 0$ near the boundary of $\Om$. We prove the existence of a weak
solution in $W^{1,p}_{loc}(\Om)$ and its behaviour near the boundary for
$\ba<p$. Consequently, we obtain optimal Sobolev regularity of weak solutions.
By establishing the comparison principle, we prove the uniqueness of weak
solution for the case $\ba<2-\frac{1}{p}$. Subsequently, for the case $\ba\ge
p$, we prove the non-existence result. Moreover, we prove H\"older regularity
of the gradient of weak solution to a more general class of quasilinear
equations involving singular nonlinearity as well as lower order terms (see
\eqref{Prb}). This result is completely new and of independent interest.
  In addition to this, we prove H\"older regularity of minimal weak solutions
of $(P)$ for the case $\beta+\delta\geq 1$ that has not been fully answered in
former contributions even for $p$-Laplace operators.
