---
title: A debiasing approach to microseismic
runninghead: Debiasing approach 
author: |
    Shashin Sharan^1^\*, Rajiv Kumar^1^, Rongrong Wang^2^, and Felix J. Herrmann^3^ \
    ^1^ Department of Earth & Atmoshpheric Sciences, Georgia Institute of Technology \
    ^2^ Michigan State University \
    ^3^ School of Computational Science and Engineering, Georgia Institute of Technology \
bibliography:
	- abstract.bib   
---


## Abstract:

Microseismic data is often used to locate fracture locations and their origin in time created by fracking. Although surface microseismic data can have large apertures and is easier to acquire than the borehole data, it often suffers from poor signal to noise ratio (S/R). Poor S/R poses a challenge in terms of estimating the correct location and source-time function of a microseismic source. In this work, we propose a denoising step in combination with a computationally cheap debiasing based approach to locate microseismic sources and to estimate their source-time functions with correct amplitude from extremely noisy data. Through numerical experiments, we demonstrate that our method can work with closely spaced microseismic sources with source-time functions of different peak amplitudes and frequencies. We have also shown the ability of our method with the smooth velocity model.

\vspace*{-0.5cm}
## Introduction
\vspace*{-0.25cm}

To make unconventional reservoirs economical for the production of oil & gas, hydraulic fracturing is a common practice adapted by the oil & gas industry. During hydraulic fracturing fractures are created, which give rise to microseismic events. To make drilling decisions and to prevent hazardous situations, we need accurate information on the location and temporal evolution of these fractures. Because microseismic waves carry important information about fracture's location and origin time, the microseismic data recorded at surface or along a monitor well is often used to locate these fractures [@maxwell2014hydraulic].

Because of the operational ease and option to cover wide aperture, surface receivers are widely used [@peter2010reservoir; @lakings2006SEGsurface]. But the microseismic data recorded along the surface comes at a cost of poor signal to noise ratio (S/R) in comparison to the data recorded along a monitor well. This is because we record more ambient noise [@arani2012SEGanalysis]. Moreover, microseismic waves suffer from attenuation while travelling large distance from subsurface to the surface receivers [@maxwell2013CSEGmicro], which makes them difficult to observe in noisy data. Thus, low S/R of surface microseismic data poses a big challenge in terms of estimation of accurate location and the origin time of microseismic sources. For example, travel-time picking based methods rely on accurate picking of first arrivals of P and S-phases. When the noise levels are high, it becomes difficult to accurately pick these first arrivals [@bolton2001travel]. Sometimes, the weak signal is not even visible to be picked. 

In @sharan2016SEGsparsity & @sharan2017SEGhigh, we proposed a computationally efficient sparsity promotion based method to invert for the microseismic source wavefield from which we extract locations and source-time functions of closely spaced microseismic sources. Our method works with noisy data with S/R as low as ``-1`` dB, but performs poorly as the S/R decreases further. Moreover, while our sparsity-promotion based method is able to give us a good estimate of the shape of the source-time function, its amplitude is often incorrect. To overcome these limitations, we propose a debiasing approach to handle the noise and correctly estimate the amplitude of the source-time function. Our proposed approach consists of two main steps. The first step involves curvelet based denoising along with sparsity promotion based microseismic source inversion to detect the location of the microseismic sources. Since noise and signal has different morphological behaviour in curvelet domain [@candes2000curvelets] it is easy to separate the noise in curvelet domain [@herrmann2008GJInps; @neelmani2008GTLECohere; @kumar2017EAGEdha]. In the second step, we use the estimated source location and perform a wave-equation based debiasing step to get the source-time function with correct amplitudes. 
     
This paper is organized as follows. We first discuss the challenges in terms of detecting the location and estimating the source-time function of microseismic sources in the presence of strong incoherent noise in the observed data. Next, we explain the basics of curvelet transform and the steps we are using to denoise the noisy observed data. Subsequently, we explain the debiasing step to get the correct amplitude of the source-time function. Finally, we show the efficacy of the proposed approach on a noisy dataset generated on a complex subset of the Marmousi model [@brougois1990marmousi]. 

\vspace*{-0.5cm}
## Methodology
\vspace*{-0.25cm}

Fracturing of rocks during fracking causes emission of microseismic waves. The microseismic events causing this emission are mostly localized along the fracture tips. Therefore, we assume these microseismic sources to be sparse in space. In @sharan2016SEGsparsity, we exploited the fact that microseismic sources are sparse in space and have finite energy along time to solve 

```math #eqBPDN
\min_{\vector{Q}} \quad \|\vector{Q}\|_{2,1} \quad \text{s.t.}\quad\|\mathcal{F}[\vector{m}] (\vector{Q}) - \vector{d}\|_2
\leq \epsilon,
```
where ``\vector{Q} \in \R^{n_x \times n_t}``, with ``n_x`` being the size of spatial grid and ``n_t`` being the number of time samples, is a matrix representing the complete microseismic source field containing the spatial temporal distribution of different microseismic sources ---i.e. the ``(i,j)`` entry in  ``\vector{Q}_{i,j} = {q}(x_i,t_j)``.  ``\mathcal{F}[\vector{m}]=\vector{P}{\mathcal{A}[\vector{m}]}^{-1}`` is the linear operator modeling the ``2``D time-domain acoustic wave-equation. The linear operator ``\vector{P}`` restricts the wavefield to the receivers. ``\mathcal{A}[\vector{m}]`` is the ``2``D finite-difference time stepping operator parametrized by the squared slowness ``\vector{m}`` of the medium. The minimizaton problem [#eqBPDN] aims to find such a ``\vector{Q}``, which has a minimum ``\ell_1``-norm in space and has a minimum ``\ell_2``-norm in time while fitting the observed data ``\vector{d}`` within the noise level ``\epsilon``. As mentioned earlier, sparseness of microseismic source wavefield ``\vector{Q}`` in space justifies the choice of ``\ell_1``-norm in space and finite energy of these microseismic sources justifies the choice of ``\ell_{2}``-norm in time in Equation [#eqBPDN]. 

Problem [#eqBPDN] has a form very similar to the classic Basis Pursuit Denoising (BPDN) problem [@chen1998atomic; @berg2008SJSCpareto]. In @sharan2016SEGsparsity, we proposed a new algorithm tailored to solve a slightly modified version of the problem [#eqBPDN] for the situations when the forward modeling operator is ill conditioned and computationally expensive. We will discuss about this new algorithm in the next section.  

### Linearized Bregman algorithm

Motivated by the recent successful application of linearized Bregman algorithm [@yin2008bregman; @dirk2014LBR] to solve sparsity promoting least-squares migration problem [@herrmann2015EAGEfast], in @sharan2016SEGsparsity, we proposed to solve

```math #eqSrcest
\min_{\vector{Q}} \quad \|\vector{Q}\|_{2,1} + \frac{1}{{2}{\mu}}\|\vector{Q}\|_F^2 \quad \text{s.t.}\quad \|\vector{\mathcal{F}}[\vector{m}] (\vector{Q}) - \vector{d}\|_2
\leq \epsilon,
```
which is strongly convex and a relaxed form of the original classic BPDN problem [#eqBPDN]. ``\|.\|_F`` is the Frobenius norm and ``\mu`` acts as a trade-off parameter between sparsity given by ``\ell_{2,1}``-norm term and the Frobenious norm term. When ``\mu \uparrow \infty``, then Equation [#eqSrcest] is equivalent to solving the original BPDN problem [#eqBPDN]. Solving Equation [#eqSrcest] can be achieved through a simple algorithm with few tuning parameters [@sharan2016SEGsparsity]. We estimate the location of microseismic sources as outliers in the intensity plot calculated as ``\vector{I}(\vector{x})=\mathrm{vec}^{-1}\left(\sum_{t}\mid\vector{Q}(\vector{x},{t})\mid\right)`` from the inverted source field ``\vector{Q}``, where ``\mathrm{vec}^{-1}(\cdot)`` reshapes a vector into its original matrix form. The temporal variation of the inverted source field ``\vector{Q}`` at the estimated source locations give the source-time function of microseismic sources.   

To avoid fitting noise in the data, every iteration of linearized Bregman algorithm involves projecting the data residual ``\vector{r}`` (i.e. difference between predicted and observed data) on the ``\ell_{2}`` ball of size ``\epsilon`` [@sharan2016SEGsparsity]

```math #project
\Pi_\epsilon (\vector{r}) = \max \left\{{0,1-\frac{\epsilon}{\|\vector{r}\|}}\right\}\,\vector{r}.
```
Linearized Bregman algorithm performs well in locating microseismic sources and estimating their source-time function when the data has low to moderate levels of noise. But this is not always the case, microseismic data can be very noisy--- i.e. ``\epsilon \gg \|\vector{r}\|``. Higher value of ``\epsilon`` implies that the projection of data residual in Equation #project will give a vector with all zeros. Therefore, linearized Bregman fails to update the source field ``\vector{Q}`` at every iteration. By using a smaller value of ``\epsilon`` instead of the actual noise level, projection in the Equation #project works and we can get source field ``\vector{Q}`` updated in every iteration of linearized Bregman algorithm. But using a smaller value of ``\epsilon`` instead of actual noise level in data means linearized Bregman algorithm will invert for such a source field ``\vector{Q}`` which will also fit the noise. Hence, the inverted source field will have many false sources as we observe this phenomenon in the numerical experiment section. To avoid the above mentioned situations with very noisy data, we propose to incorporate curvelet based denoising step prior to applying linearized Bregman algorithm to invert for the location followed by a wave-equation based debiasing approach to get the source-time function with correct amplitude.   

### Curvelet based denoising

Microseismic signals recorded by the surface receivers are very weak in amplitude and are often contaminated with ambient noise that have similar or higher amplitude level than the amplitude of the microseismic signal present in the data [@arani2012SEGanalysis]. Also, the frequency range of the ambient noise is very similar to that of the microseismic signal [@onge2011noise]. This makes signal and noise separation very difficult and eventually causes problems in detecting the microseismic sources and estimating their source-time functions. Curvelet transform is a multi-scale and multi-directional transform [@candes2000curvelets], that maps seismic data into angular wedges of different scales in the ``2``D Fourier domain (Figure [#curvelets]). This property of the curvelet transform helps in separating signals components based on their location, dip and scaling in the transform domain. Therefore, curvelet transform has been succefully used for incoherent [@herrmann2008GJInps; @neelmani2008GTLECohere] and coherent noise attenuation [@kumar2017EAGEdha; @lin2013GEOPrepsi].     

#### Figure: {#curvelets}
![](Figure/curvelet.png){#acq1a width=50%}
![](Figure/curvelet1.png){#acq1b width=50%}

: Properties of the curvelet transform (Source: @herrmann2008GJInps)

Motivated by prior successful application of curvelets, we propose following steps for denoising:

### Algorithm: Curvelet based denoising {#alg:denoise}
| 1. Noisy data ``\vector{d}``, forward curvelet transform operator ``\vector{C}``, Threshold parameter ``\lambda``  `//Input`
| 2. ``\vector{b} = \vector{C}\vector{d}`` `//Forward curvelet transform`
| 3. ``[\vector{sb},\vector{idx}] = \textrm{Sort}(\mid\vector{b}\mid)`` `//Sorting in descending order`
|      where ``\vector{idx}`` stores the indices of sorted curvelet coefficients
| 4. ``(\vector{e})_{h} = \sqrt{\frac{\sum_{i=1}^{h}\vector{sb}_{i}^2}{\sum_{i=1}^{l}\vector{sb}_{i}^2}}`` `//normalized cumulative energy`
|      where ``{l}`` is the length of ``\vector{sb}``
| 5. Find the smallest index ``{p}`` such that ``(\vector{e})_{p} \geq \lambda``
| 6. ``\vector{R} = \vector{C}^\top(\vector{idx}(1:{p}),:)`` `//New inverse curvelet transform operator`
| 7. ``\vector{b}_{dn} = (\vector{R}^\top\vector{R})^{-1}\vector{R}^\top\vector{d}`` `//Solving the normal equation`
| 8. ``\vector{d}_{dn} = \Re(\vector{R}\vector{b}_{dn})`` `//denoising`
Caption: Denoising with curvelets.

**Line 2** in the Algorithm #alg:denoise corresponds to forward curvelet transform of the noisy microseismic data ``\vector{d}`` (Figure [#demoa]) in the curvelet domain (Figure [#demob]). The indices in Figure [#demob] are arranged from coarse to fine scale. **Line 3** corresponds to sorting of the absolute value of curvelet coefficients ``\vector{sb}`` in descending order. We also store the indices of the sorted curvelet coefficients in ``\vector{idx}``. **Line 4** corresponds to computing the square root of normalized cumulative energy of the sorted curvelet coefficients. **Line 5** corresponds to finding the smallest index ``{p}`` in vector ``\vector{e}`` at which the square root of the normalized cumulative energy of sorted curvelet coefficient exceeds or is equal to the threshold ``\lambda``. In **line 6**, we form a subset ``\vector{R} \subseteq \vector{C}^\top`` of the inverse curvelet transform operator. Columns of ``\vector{R}`` correspond to the curvelet coefficients whose square root of the normalized cumulative energy in ``\vector{e}`` is greater than or equal to the threshold ``\lambda``. **Line 7** corresponds to solving the normal equation to get the debiased and denoised curvelet coefficients ``\vector{b}_{dn}`` effectuated by the new inverse curvelet transform operator ``\vector{R}``. Debiasing neutralizes the shrinkage effect of thresholding and preserves energy. Figure [#democ] shows absolute value of the denoised and debiased curvelet coefficients ``\vector{b}_{dn}`` mapped to the corresponding location of the noisy curvelet coefficients ``\vector{b}``. **Line 8** corresponds to taking inverse curvelet transform of the denoised curvelet coefficients ``\vector{b}_{dn}`` and taking its real part to get the denoised microseismic data (Figure [#demod]). We choose the threshold parameter as large as possible but for which we do not see any primary leakage in the difference plot.     

The curvelet based denoising involves very few forward and inverse curvelet transform, which makes the proposed denoising method computationally cheap. To detect the location of microseismic sources from the denoised data ``\vector{d}_{dn}`` in a computationally efficient manner, we use the accelerated version of linearized Bregman algorithm [@sharan2017SEGhigh] along with a ``2``D left preconditioner [@herrmann2009GEOPcbm]. 

#### Figure: {#demo}
![](Figure/Shot_high_noise_Marm.png){#demoa width=50%}
![](Figure/noisy_coefficient.png){#demob width=50%}\
![](Figure/denoised_coefficient.png){#democ width=50%}
![](Figure/Shot_h_denoise_Marm.png){#demod width=50%} 

: Curvelet denoising schematic. *(a)* Noisy microseismic data with S/R = ``-5.70`` dB. *(b)* Absolute value noisy curvelet coefficient. *(c)* Absolute value debiased curvelet coefficients. *(d)* Denoised microseismic data with S/R = ``5.3`` dB.  

### Debiasing of the source-time function

Given the location of microseismic sources, next step is to estimate the correct amplitude of source-time function. To achieve this, we use the forward modeling operator ``\vector{\mathcal{F}}[\vector{m}]`` and estimate source locations to fit the noisy data ``\vector{d}`` within some tolerance level. We use noisy data to avoid any kind of amplitude errors introduced in the approximated data by denoising. We now solve a debiasing problem using least squares as

```math #eq:matee
\tilde{\vector{W}} = \argmin_{\vector{W}\in\R^{n_t \times n}} \quad \|\mathcal{F}[\vector{m}](\vector{H}\vector{W}^\top) - \vector{d}\|,
```

where ``\vector{H}\in\R^{n_x \times n}`` is a matrix, with ``{n}`` being the number of detected microseismic sources, whose ``{i}^{th}`` column is ``\vector{h}_{i}``,which corresponds to the location of the ``{i}^{th}`` source. ``\vector{h}_{i}`` is a spatial delta function ``\delta(\vector{x}-\vector{x}_{i})`` with ``\vector{x}_{i}`` being the location of ``i^{th}`` microseismic source. Equation [#eq:matee] solves for the unknown matrix ``\vector{W}`` whose ``{i}^{th}`` column corresponds to the source-time function of ``{i}^{th}`` microseismic source. We run only a few iterations of the unconstrained problem [#eq:matee] to avoid overfitting the noise in the data.

\vspace*{-0.5cm}
## Numerical Experiments
\vspace*{-0.25cm}


To demonstrate the effectiveness of our method for data with high noise level (S/R = ``-7.30`` dB), we performed a numerical experiment. We used ``2``D acoustic finite-difference modeling code [@louboutin2017fwi] to generate microseismic data of record length ``1.0\, \mathrm{s}``. To make the experimental setup more realistic, we used ``5`` microseismic sources of different amplitudes (differ by a factor of 2) and dominant frequencies ``(30.0\,``, and  ``25.0)\, \mathrm{Hz}`` activating at a small time interval with overlapping source-time functions to generate the microseismic data. Because of its geological complexity, we chose a part of Marmousi model with dimensions ``3.15\,\mathrm{km} \times 1.08\,\mathrm{km}`` (``631 \times 217`` points) to perform the experiment. We place ``5`` microseismic sources (indicated by black dots in Figure [#acqa]) in low velocity layer to generate the data. The adjacent sources are separated by half a dominant wavelength. We use surface receivers placed at a depth of ``20.0\, \mathrm{m}`` from the top surface to record the data. To get noisy data (Figure [#estsrca]) we add random noise (``5.0\, \mathrm{Hz}`` to ``40.0\, \mathrm{Hz}``) to the noise free data. We use kinematically correct smooth velocity model to invert for microseismic source field in the experiment. As expected, our method performs poorly without curvelet denoising and gives an intensity plot (Figure [#acqb]) that is not very informative. This is because of the presence of lots of false sources in the estimated intensity plot (Figure [#acqb]). Therefore, we apply the proposed curvelet based denoising steps to the noisy data (Figure [#estsrca]) to get denoised data (Figure [#estsrcb]) with improved S/R of ``3.5`` dB. The difference plot (Figure [#estsrcc]) between noisy (Figure [#estsrca]) and the denoised data (Figure [#estsrcb]) shows that we do not loose any coherent signal with the proposed denoising method. With only ``10`` iterations of accelerated version of linearized Bregman algorithm, we are now able to locate all the ``5`` microseismic sources (Figure [#estsrcd]) from this denoised data (Figure [#estsrcb]). The white colour crosses in the estimated intensity plot correspond to the actual location of microseismic sources. All the outliers are located near the actual location of microseismic sources. To get the correct source-time function (blue colour plot in Figures [#estwaveleta] to [#estwavelete]), we perform debiasing by least-squares using the estimated source location. We use the noisy data to perform this debiasing. Denoising helps us to get the correct source location and the debiasing step helps us to recover the source-time function with correct scaling. We further compare the source-time function (blue colour plot in Figures [#estwaveleta] to [#estwavelete]) estimated by proposed approach to the source-time function estimated by @sharan2017SEGhigh (red colour plot in Figures [#estwaveleta] to [#estwavelete]). For visualization purpose, we scale the wavelets displayed in red color by a factor of ``40``. Thus, the proposed method can estimate the location and source-time function with correct scaling for microseismic data acquired in extremely noisy environment.          

### Figure: {#acq}

![](Figure/True_Vel_Marm_box.png){#acqa width=50%}
![](Figure/Est_src_l21_mu9emin4_10iter_very_high_noise_P1_no_denois_epsp1.png){#acqb width=50%}

: *(a)* Acquisition geometry with velocity model. Inverted yellow colour triangles indicate receivers buried at a depth of ``20.0\, \mathrm{m}`` & separated by ``10.0\, \mathrm{m}``. Black dots indicate the location of two microseismic sources. *(b)* Estimated intensity plot from noisy data without denoising. White colour dots indicate actual location of microseismic sources.


### Figure: {#estsrc}

![](Figure/Shot_very_high_noise_Marm.png){#estsrca width=50%}
![](Figure/Shot_vh_denoise_Marm.png){#estsrcb width=50%}\
![](Figure/Shot_vh_diff_Marm.png){#estsrcc width=50%}
![](Figure/Est_src_l21_mu9emin4_10iter_very_high_noise_P1_denois.png){#estsrcd width=50%}

: Noisy microseismic data and estimated intensity plots(zoomed): noisy data with *(a)* S/R = ``-7.3`` dB. *(b)* Denoised data using curvelet based denoising with improved S/R of ``3.5`` dB. *(c)* Data difference plots after denoising. *(d)* Estimated intensity plot. White colour crosses indicate the true location of microseismic sources.


### Figure: {#estwavelet}

![](Figure/Wavelet_comparison_lsqr_loc1.png){#estwaveleta width=50%}
![](Figure/Wavelet_comparison_lsqr_loc2.png){#estwaveletb width=50%}\
![](Figure/Wavelet_comparison_lsqr_loc3.png){#estwaveletc width=50%}
![](Figure/Wavelet_comparison_lsqr_loc4.png){#estwaveletd width=50%}\
![](Figure/Wavelet_comparison_lsqr_loc5.png){#estwavelete width=50%}

: Source-time function comparison: Comparison of the true source-time functions (solid green) with source-time function (blue color) estimated by proposed method. We also perform comparison with the source-time function (amplified by ``40`` times) estimated using the approach proposed in @sharan2017SEGhigh (solid magenta) at locations *(a)* ``1``, *(b)* ``2``, *(c)* ``3``, *(d)* ``4``, *(e)* ``5`` from LtoR in Figure [#acqa]. Dominant frequency of source-time functions at (from LtoR) locations ``1`` and ``2`` is ``25.0\, \mathrm{Hz}``, at locations ``3``, ``4`` and ``5`` dominant frequency is ``30.0\, \mathrm{Hz}``.


\vspace*{-0.5cm}
## Conclusions
\vspace*{-0.25cm}

We proposed a debiasing based approach to estimate the location and source-time function of microseismic sources with correct amplitude from data with very low S/R. We showed ability of the proposed method to resolve microseismic events even when the sources are spatially close and have overlapping source-time functions. The proposed method is also computationally cheap as it requires very few forward and backward curvelet transforms along with few iterations of the accelerated version of linearized Bregman. Also, the source-time function estimation requires only a very few least-squares iterations. In future, we would like to apply PCA based denoising techniques to deal with different types of noise such as ground roll, source side noise etc.
  

\vspace*{-0.5cm}
## Acknowledgement
\vspace*{-0.25cm}

This research was carried out as part of the SINBAD II project with the support of the member organizations of the SINBAD Consortium.




```math_def
\def\argmin{\mathop{\rm arg\,min}}
\def\vec{\mbox{$\mathrm{vec}$}}
\def\ivec{\mbox{$\mathrm{vec}^{-1}$}}
%\newcommand{\rr}{{{r}}}
\newcommand{\m}{{\mathsf{m}}}
\newcommand{\PsDO}{\mbox{PsDO\,}}
\newcommand{\Id}{\mbox{$\tensor{I}\,$}}
\newcommand{\R}{\mbox{$\mathbb{R}$}}
\renewcommand{\C}{\mbox{$\mathbb{C}$}}
\newcommand{\E}{\mbox{$\mathbb{E}$}}
\newcommand{\Z}{\mbox{$\mathbb{Z}$}}
\newcommand{\DE}{:=}
\newcommand{\Order}{\mbox{${\cal O}$}} \def\bindex#1{{\mathcal{#1}}}
%\newcommand{\Id}{\mbox{$\tensor{I}\vector{d}\,$}}
\def\pector#1{\mathrm{\mathbf{#1}}}
\def\cector#1{#1}
\def\censor#1{#1}
\def\vector#1{\mathbf{#1}}
\def\fvector#1{{\widehat{\vector{#1}}}}
\def\evector#1{{\widetilde{\vector{#1}}}}
\def\pvector#1{{\breve{\vector{#1}}}}
\def\pector#1{\mathrm{#1}}
\def\ctensor#1{\bm{\mathcal{#1}}}
%\def\tensorm#1{{#1}}
\def\tensorm#1{\bm{#1}}
%\def\tensor#1{\bm{#1}}7
\def\tensor#1{\vector{#1}}
%\def\cector#1{{\vector{\underline{#1}}}}
\def\hensor#1{\tensor{#1}}
\def\uensor#1{\underline{\bm{#1}}}
%\def\hensor#1{{\boldsymbol{\mathsf{#1}}}}
\def\hector#1{\vector{#1}}
%\def\hector#1{{\boldsymbol{\mathsf{#1}}}}
% \def\hfector#1{\hat{\boldsymbol{\mathsf{#1}}}}
% \def\ctensor#1{{\tensor{\underline{\underline{\mathsf{#1}}}}}}
\DeclareMathOperator\minim{\hbox{minimize}}
\def\ftensor#1{{\widehat{\tensor{#1}}}}
\def\calsor#1{{\boldsymbol{\mathcal{#1}}}}
\def\optensor#1{{\boldsymbol{\mathcal{#1}}}}
\def\hvector#1{\hat{\boldsymbol{\mathbf{#1}}}}
\def\uvector#1{{\underline{\vector{#1}}}}
\def\utensor#1{{\underline{\tensor{#1}}}}
%\
\newcommand{\dd}{\ensuremath{\mathrm{d}}}
\newcommand{\blockstack}{\operatorname{blockstack}}
\newcommand{\diag}{\operatorname{diag}}
\newcommand{\rs}[1]{\mathstrut\mbox{\scriptsize\rm #1}}
\newcommand{\rr}[1]{\mbox{\rm #1}}
\usepackage{algorithm2e}
\def\minim{\mathop{\hbox{minimize}}}
```