The constraints can be rewritten by replacing the squares with absolute values, then you have a program like

$\operatorname{argmin}_{\mathbf{X}} \|\mathbf{D}\mathbf{X}\|_F^2 + \|\mathbf{X} - \mathbf{Z}\|_F^2$

$\text{s.t.}\quad \forall i,j \quad |\mathbf{d}_{01} \mathbf{X}_i| = |\mathbf{d}_{01} \mathbf{X}_j|, \text{ and } |\mathbf{d}_{34} \mathbf{X}_i| = |\mathbf{d}_{34} \mathbf{X}_j| $

where the $\mathbf{X}_i$ constitute the columns of $\mathbf{X},$ the matrix $\mathbf{D}$ takes the differences of the consecutive columns of $\mathbf{X},$ and $\mathbf{d}_{01}, \mathbf{d}_{34}$ are appropriately defined row vectors.

Then you can use the standard trick of replacing $|x|$ with $x_{+} + x_{-}$ subject to $x = x_{+} - x_{-}$ and $x_{+}, x_{-} \geq 0$ to get the equivalent program

$\operatorname{argmin}_{\mathbf{X}, \mathbf{e}^{+}, \mathbf{e}^{-},\mathbf{f}^{+}, \mathbf{f}^{-}} \|\mathbf{D}\mathbf{X}\|_F^2 + \|\mathbf{X} - \mathbf{Z}\|_F^2$

$\text{s.t.}\quad \forall i: \quad \mathbf{d}_{01} \mathbf{X}_i = e^{+}_i - e^{-}_i, \text{ and } \mathbf{d}_{34} \mathbf{X}_i = f^{+}_i - f^{-}_i,$

$\quad \quad \forall i \neq j: e^{+}_i + e^{-}_i = e^{+}_j + e^{-}_j, \text{ and } f^{+}_i + f^{-}_i = f^{+}_j + f^{-}_j,$

$\quad \quad \mathbf{e}^{+}, \mathbf{e}^{-}, \mathbf{f}^{+}, \mathbf{f}^{-} \succeq 0$

By considering the Lagrangian, you can see that for each i, at least one of the $e^{\pm}_i$ will be zero, and likewise for the $f^{\pm}_i,$ so this is indeed equivalent to the original program.

By vectorizing $\mathbf{X}$, you can write the objective as a convex quadratic in a vector $\mathbf{x}$, and a little massaging of the constraints gives you a convex quadratic program in standard form.

In terms of actually solving the final QP, I'm not sure what the state of art is, especially performance usually depends on the properties of your problem (sparsity or structure in the hessian of your quadratic, whether you have only equality of inequality constraints, etc.). You'd need a subject matter expert to get a good answer on that one. I'd personally try TFOCS if you have access to Matlab, since it's a first order solver and looks relatively easy to use (I have not used it myself).

The constraints can be rewritten by replacing the squares with absolute values, then you have a program like

$\operatorname{argmin}_{\mathbf{X}} \|\mathbf{D}\mathbf{X}\|_F^2 + \|\mathbf{X} - \mathbf{Z}\|_F^2$

$\text{s.t.}\quad \forall i,j \quad |\mathbf{d}_{01} \mathbf{X}_i| = |\mathbf{d}_{01} \mathbf{X}_j|, \text{ and } |\mathbf{d}_{34} \mathbf{X}_i| = |\mathbf{d}_{34} \mathbf{X}_j| $

where the $\mathbf{X}_i$ constitute the columns of $\mathbf{X},$ the matrix $\mathbf{D}$ takes the differences of the consecutive columns of $\mathbf{X},$ and $\mathbf{d}_{01}, \mathbf{d}_{34}$ are appropriately defined row vectors.

Then you can use the standard trick of replacing $|x|$ with $x_{+} + x_{-}$ subject to $x = x_{+} - x_{-}$ and $x_{+}, x_{-} \geq 0$ to get the equivalent program

$\operatorname{argmin}_{\mathbf{X}, \mathbf{e}^{+}, \mathbf{e}^{-},\mathbf{f}^{+}, \mathbf{f}^{-}} \|\mathbf{D}\mathbf{X}\|_F^2 + \|\mathbf{X} - \mathbf{Z}\|_F^2$

$\text{s.t.}\quad \forall i: \quad \mathbf{d}_{01} \mathbf{X}_i = e^{+}_i - e^{-}_i, \text{ and } \mathbf{d}_{34} \mathbf{X}_i = f^{+}_i - f^{-}_i,$

$\quad \quad \forall i \neq j: e^{+}_i + e^{-}_i = e^{+}_j + e^{-}_j, \text{ and } f^{+}_i + f^{-}_i = f^{+}_j + f^{-}_j,$

$\quad \quad \mathbf{e}^{+}, \mathbf{e}^{-}, \mathbf{f}^{+}, \mathbf{f}^{-} \succeq 0$

By considering the Lagrangian, you can see that for each i, at least one of the $e^{\pm}_i$ will be zero, and likewise for the $f^{\pm}_i,$ so this is indeed equivalent to the original program. Concisely: at a presumed optimal point you can always increase over the value $L(\mathbf{X}^\star, \mathbf{e}^{+,\star}, \mathbf{e}^{-,\star}, \mathbf{f}^{+,\star}, \mathbf{f}^{-,\star}, \boldsymbol{\lambda}, \boldsymbol{\nu})$ by moving all the mass onto exactly one of $\mathbf{e}^{\pm}$ and $\mathbf{f}^{\pm}.$

By vectorizing $\mathbf{X}$, you can write the objective as a convex quadratic in a vector $\mathbf{x}$, and a little massaging of the constraints gives you a convex quadratic program in standard form.

In terms of actually solving the final QP, I'm not sure what the state of art is, especially since performance usually depends on the properties of your problem (sparsity or structure in the Hessian of your quadratic, whether you have only equality or inequality constraints, etc.). You'd need a subject matter expert to get a definitive answer on that one. I'd personally try TFOCS if you have access to Matlab, since it's a first order solver and looks relatively easy to use (I have not used it myself).