--- title: High resolution fast microseismic source collocation and source time function estimation runninghead: Sparsity-promoting source estimation author: | Shashin Sharan^1^\*, Rongrong Wang^1^, and Felix J. Herrmann^1^ \ ^1^ Seismic Laboratory for Imaging and Modelling (SLIM), University of British Columbia bibliography: - abstract.bib --- ## Abstract: Sparsity promotion based joint microseismic source collocation and source time function estimation, using Linearized Bregman algorithm, although simple to implement, suffers from slow convergence. This is due to the fact that Linearized Bregman algorithm has only first order of convergence. In this work, we propose to accelerate the existing Linearized Bregman algorithm using the L-BFGS algorithm. Without any initial guess for the source location or source time function, our method is able to estimate the source location and source time function for kinematically correct smooth velocity model. We demonstrate the effectiveness of our method for multiple sources spaced within half a wavelength. We also compare our results with Linearized Bregman based method in ``2.5``D instead of ``2``D. \vspace*{-0.5cm} ## Introduction \vspace*{-0.25cm} There has been increased attention towards unconventional reservoirs in last couple of decades. Unconventional reservoirs require some kind of stimulation in the form of hydraulic fracturing to make them economically viable for production by connecting weakly connected oil & gas traps through the fractures. For production purposes and to prevent any kind of hazardous situations (e.g. interference of these fractures with pre-existing wells and faults) it is extremely important to estimate the fracture locations accurately together with other attributes such as source time function, stress orientation along these fractures, etc. Among all the available methods, microseismic based methods provide most accurate way to estimate the above mentioned attributes of fractures [@maxwell2014hydraulic]. As fractures caused by hydraulic fracturing propagate through the subsurface rocks, they generate microseismic waves that are recorded by receivers along the surface or receivers along monitor wells. The recorded microseismic data contains information about the source attributes such as source location, source time function, source mechanism etc. One exploits this microseismic data to find out different source attributes, which is extremely helpful for drilling decisions, hazard prevention etc. Traditional methods based on travel time inversion require travel time picking and identification of phases such as P and S wave components [@thurber2000earth; @waldhauser2000double]. These methods estimate the location but not the source time function. Recently developed methods based on imaging [@george1982det; @dirk2005reverse; @sun2015SEGitp; @nakata2016reverse] do not require travel time picking but only estimate source location, again without the source time function. @kad2015SEGmee suggested a FWI based method to estimate both source location and source time function in an alternating fashion. @wang2016SEGmicro proposed to use FWI to invert for source location, source time function and the velocity model in an alternating fashion. Both of these works assumed that the spatial and temporal component of a source is separable. Moreover, they showed results only for a single source. However, in real life there can be more than one microseismic sources. @sharan2016SEGsparsity showed that in order to use FWI with the assumption of separable source location and source time function, one need to know the number of sources beforehand. For hydraulic fracturing scenario the knowledge of number of microseismic sources is not a trivial problem. Also, FWI based method suggested by @kad2015SEGmee performs poorly when both the source location and source time function are unknown, which is mostly the situation for real life scenario. In @sharan2016SEGsparsity, the authors exploited the fact that wave equation itself is an excellent sparsity promoting transform to simultaneously estimate the microseismic source location & source time function without any prior assumptions about the number of sources or any initial guess about the source location and source time function. They used the Linearized Bregman algorithm to jointly estimate the microseismic source location & source time function. Although, this algorithm is simple to implement, it suffers from slow convergence because of the 1st order convergence. As each iteration of the algorithm required two PDE solves (Forward and adjoint wave-equation solve in time domain), using the Linearized Bregman algorithm can be too expensive for a complex problem. @huang2013accel and @yin2010analysis showed that the Linearized Bregman algorithm is equivalent to solving a gradient descent of the dual problem and therefore they proposed to use a Limited memory BFGS [@liu1989on] kind of method to accelerate the convergence of original Linearized Bregman algorithm from ``1^{st}`` order to ``2^{nd}`` order of convergence. Motivated by this idea, in this work, we propose to use Linearized Bregman algorithm with LBFGS to attain faster convergence for the joint microseismic source-collocation and source time function estimation problem. We also extend our proposed method to the situations involving 2.5 D wave propagation instead of 2D wave propagation. We then demonstrate the effectiveness of our method for complex geology and closely spaced sources with different source-time functions with different dominant frequencies. \vspace*{-0.5cm} ## Methodology \vspace*{-0.25cm} Given a sufficiently accurate velocity model, @sharan2016SEGsparsity exploited the fact that the wave equation acts as a sparsity promoting "transform" by focusing wavefields on their correct source locations. Assuming sources are localized in space and have finite energy along time, @sharan2016SEGsparsity proposed to solve ```math #eqMoti \min_{\vector{U}} \quad \|\mathcal{A}[\vector{m}]\ (\vector{U})\|_{2,1} \quad \text{s.t.} \quad \|\mathcal{P}(\vector{U}) - \vector{d}\|_2 \leq \epsilon. ``` Equation [#eqMoti] aims to minimize the ``\|.\|_{2,1}`` norm of action of finite difference time stepping operator ``\mathcal{A}[\vector{m}]``, parametrized by squared slowness ``\vector{m}``, on spatio-temporal wavefield ``\vector{U}`` under the condition that receiver restriction operator ``\mathcal{P}`` restricts the wavefield ``\vector{U}`` to fit the observed microseismic data ``\vector{d}`` within ``\epsilon``, the two-norm of noise level. The ``\|.\|_{2,1}`` norm was used to exploit the fact that sources have finite energy along time but are sparse along space. The operators ``\mathcal{A}[\vector{m}] (\vector{U})`` and ``\mathcal{P}(\vector{U})`` are linear in terms of wavefield ``\vector{U}``. @sharan2016SEGsparsity proposed to make the above optimization problem more tractable by a simple change of variable ``\vector{Q}=\mathcal{A}[\vector{m}] (\vector{U})``, where ``\vector{Q}`` is a matrix containing spatial temporal distribution of all the sources---i.e. the ``(i,j)`` entry in ``\vector{Q}_{i,j} = {q}(x_i,t_j)``. So equation [#eqMoti] becomes: ```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 ``\mathcal{F}[\vector{m}]=\mathcal{P}{\mathcal{A}[\vector{m}]}^{-1}`` is the forward modeling operator which is a linear operator in terms of ``\vector{Q}`` . Equation [#eqBPDN] is similar to the classic basis pursuit denoising (BPDN) problem. Motivated by recent successful application of the linearized Bregman [@yin2008bregman; @dirk2014LBR] to least-squares migration [@herrmann2015EAGEfast], @sharan2016SEGsparsity proposed to solve slight convex relaxation of the above optimization problem. They 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, ``` where ``\|.\|_F`` is the Frobenius norm and ``\mu`` acts as a trade off parameter between the sparsity and regular two norm terms. When ``\mu \uparrow \infty``, equation [#eqSrcest] takes the form of the BPDN problem in equation [#eqBPDN]. Although Linearized Bregman is a simple ``3`` step algorithm, its iterations require solving the wave equation and its adjoint to jointly collocate the source locations and estimate the associated source signatures [@sharan2016SEGsparsity]. The cost of each iteration of Linearized Bregman depends on the model size, model dimension (``2``D or ``3``D), and on the number of time samples. For closely spaced sources within half a wavelength, we need a higher value of ``\mu`` to ensure sparse solutions. However, a higher value of ``\mu`` makes the convergence of the Linearized Bregman algorithm slow. Therefore, it will require too many iterations of Linearized Bregman, which is computationally expensive. @huang2013accel and @yin2010analysis proposed a solution to overcome the problem of slow convergence of Linearized Bregman. They proved that the Linearized Bregman method applied to an optimization problem is equivalent to solving a gradient descent of its Lagrangian dual. @huang2013accel also suggested that gradient descent accelerating methods such as Limited memory BFGS [@liu1989on] can be incorporated to solve the dual problem of the original optimization problem. By doing so, they upgraded the algorithm from ``1^{st}`` order to ``2^{nd}`` order. From their observations, we propose an algorithm to solve our ``\ell_{21}``-norm minimization problem with a ball constraint, in less number of iterations. In the original Linearized Bregman algorithm [@sharan2016SEGsparsity] we only considered the primal variable, which was the whole source field ``\vector{Q}``. Now we also work with the dual variable, which has the dimensions of the observed data, to have a better approximation of the inverse of the Hessian in the L-BFGS algorithm. The advantage of this dual formulation is that we can store the dual variable updates, because of its low memory requirement, to build the approximation of the inverse of the Hessian. Hence, by using L-BFGS and the dual variable we can accelerate the convergence without incurring any extra cost in terms of memory storage. To formulate the dual problem, we first write Equation [#eqSrcest] into the following unconstrained form: ```math #eqequiv \min_{\vector{Q}} \quad \|\vector{Q}\|_{2,1} + \frac{1}{{2}{\mu}}\|\vector{Q}\|_F^2 + \iota_{\|\vector{\mathcal{F}}[\vector{m}] (\vector{Q}) - \vector{d}\|_2\leq\epsilon}\,, ``` where ``\iota_{C}`` is the support function for the set ``C``, defined as ```math #eqiota \iota_{C}(x)=\left\{\begin{matrix}0, & x \in C, \\ \infty, & x\notin C. \end{matrix}\right. ``` Applying the fact that ``\iota_{\|\vector{\mathcal{F}}[\vector{m}] (\vector{Q}) - \vector{d}\|_2\leq\epsilon} =\sup\limits_\vector{y} \langle \vector{\mathcal{F}}[\vector{m}] (\vector{Q}) - \vector{d}, \vector{y} \rangle - \epsilon\|\vector{y}\|_2 `` to ([#eqequiv]), we obtain the following saddle point formulation ```math #eqlag \min_{\vector{Q}} \quad \|\vector{Q}\|_{2,1} + \frac{1}{{2}{\mu}}\|\vector{Q}\|_F^2 - \max_{\vector{y}} \quad \big\{\vector{y}^T(\vector{\mathcal{F}}[\vector{m}] (\vector{Q}) - \vector{d})-\epsilon\|\vector{y}\|_{2}\big\} ``` To get the dual formulation of the original problem, we can rewrite equation [#eqlag] as ```math #eqlageq \max_{\vector{y}} \quad \min_{\vector{Q}} \quad \big\{ \|\vector{Q}\|_{2,1} + \frac{1}{{2}{\mu}}\|\vector{Q}\|_F^2 - \vector{y}^T(\vector{\mathcal{F}}[\vector{m}] (\vector{Q}) - \vector{d})\big\}-\epsilon\|\vector{y}\|_{2} ``` where ``\vector{y}`` is the dual variable which has the same smaller dimension as the data ``\vector{d}``. The function inside the curly bracket is a function of ``\vector{y}``. Denoting it by ``G(\vector{y})``, the above expression ([#eqlageq]) becomes ```math #eqdual \max_{\vector{y}} \quad G(\vector{y})-\epsilon \|\vector{y}\|_2. ``` @huang2013accel explained that both the function ``G(\vector{y})`` and its derivative have close form representations. As a result, the dual problem ([#eqdual]) is unconstrained, has differentiable objective function, and has a close-form derivative. For convenience, we now rewrite the optimization problem ([#eqdual]) as ```math #dual \min_{\vector{y}} \quad -\{G(\vector{y})-\epsilon \|\vector{y}\|_2\}. ``` We solve the optimization problem ([#dual]) with the L-BFGS iteration (see algorithm [#alg:LBFGS] below). ### Algorithm: LBFGS algorithm {#alg:LBFGS} | 1. **for k=1,2,…** | 2. ``\vector{V}_{k+1}=\mathcal{F}[\vector{m}]^T(\vector{y}_{k})`` `//adjoint solve` | 3. ``\vector{Q}_{k+1}=\mathrm{Prox}_{\mu}\|.\|_{2,1}(\mu\vector{V}_{k+1})`` `//sparsity promotion` | 4. ``z_{k+1}=\|\vector{Q}_{k+1}\|_{2,1} + \frac{1}{2\mu}\|\vector{Q}_{k+1}-\mu\vector{V}_{k+1}\|_F^2`` | 5. ``\vector{y}_{k+1}=\argmin_\vector{y} -\vector{d}^T\vector{y}+\mu\|\vector{V}_{k+1}\|_F^2+\epsilon\|\vector{y}\|_2-z_{k+1}`` | 6. **end** Caption: LBFGS algorithm In algorithm [#alg:LBFGS] the operator ``\mathrm{Prox}_{\mu}\|.\|_{2,1}(\vector{C}):= \argmin_\vector{B} \quad \|\vector{C}\|_{2,1} + \frac{1}{2\mu}\|\vector{C}-\vector{B}\|_F^2`` is the proximal mapping [@combettes2011proximal] of the ``\ell_{21}`` norm. Line ``2`` in the algorithm [#alg:LBFGS] corresponds to back propagation of the current estimate of the dual variable ``\vector{y}`` in time and space to update the variable ``\vector{V}``. In the next step, we scale ``\vector{V}`` by parameter ``\mu`` followed by a "focusing" via the proximal operator to get the current estimate for the source wavefield ``\vector{Q}`` both in space and in time. Step ``4`` corresponds to updating the auxiliary variable ``z`` and step ``5`` corresponds to updating the dual variable ``\vector{y}`` as a minimizer of a problem equivalent to problem ([#dual]). Hence, we get the estimate for the source wavefield ``\vector{Q}`` as a byproduct of solving the minimization problem [#dual]. Since, the dual variable ``\vector{y}`` has the dimension of data, we can store its updates to better approximate the inverse of the Hessian for the L-BFGS algorithm. Once we solve for the source ``\vector{Q}``, we generate an intensity plot ``{I}(\vector{x})=\sum_{t}\mid\vector{Q}(\vector{x},{t})\mid`` by summing the absolute values of ``\vector{Q}(\vector{x},{t})`` along time at each grid point. The outliers in the intensity plot ``{I}(\vector{x})`` correspond to the estimated source locations, the temporal variations of the estimated source wavefield ``\vector{Q}(\vector{x},{t})`` at the estimated source locations correspond to the estimated source-time function. Figure [#conv] shows the comparison of data error variation with iterations of L-BFGS and Linearized Bregman (LBR) for a given microseismic dataset and for a given value of ``\mu``. We clearly observe that the data error decays faster for LBFGS based method in comparison to Linearized Bregman based method. Here, we choose a higher value of ``\mu`` to get sparse solutions. A higher value of ``\mu`` does not allow the source wavefields ``\vector{Q}`` to be updated in the early iterations. Therefore, the Linearized Bregman algorithm goes through initial phase of stagnation. In Figure [#conv] we also demonstrate the advantage of storing dual variable updates for improved L-BFGS convergence. We achieve improved convergence without incurring any additional memory cost, since our dual variable lives in data domain. We compare the convergence for update history size of ``{5}``, ``{20}`` and ``{40}``. We observe better convergence with increasing L-BFGS update history size. ### Figure: {#conv} ![](Figure/LBR_vs_LBFGS_conv_with_legend.png){#conv width=50%} : Data error ``\ell_{2}`` norm decay with Linearized Bregman (solid purple curve) and L-BFGS with L-BFGS history size of ``{5}`` (solid blue curve), ``{20}`` (solid red curve) and ``{40}`` (solid orange curve) \vspace*{-0.5cm} ## Numerical Experiments \vspace*{-0.25cm} We performed two different experiments to compare the efficacy of the proposed method with respect to Linearized Bregman based method [@sharan2016SEGsparsity]. For both of these experiments, we assume that microseismic sources are monopole point sources. To demonstrate the effectiveness of our method in a complex geological situation, we use a part of Marmousi model with dimensions ``3.15\,\mathrm{km} \times 1.08\,\mathrm{km}`` (``631 \times 217`` points) (Figure [#acq1a]). We use ``2.5``D time-stepping acoustic finite difference modeling code to generate microseismic data. We set surface receivers buried at depth of ``20\,\mathrm{m}`` with separation of ``10\,\mathrm{m}`` to record P-wave data. We use noise free data to simultaneously estimate source location and source time function. For both experiments we consider two closely spaced sources within half a wavelength. The first source is located at (``1.625\, \mathrm{km}`` , ``0.73\, \mathrm{km}``) and the second source is located at (``1.665\, \mathrm{km}`` , ``0.72\, \mathrm{km}``) (Figure [#acq1a]). Both sources are in low velocity zone and separated within half a wavelength. Both sources are activated at different times and with different dominant frequencies of ``30.0\, \mathrm{Hz}`` and ``25.0\, \mathrm{Hz}`` respectively. We use ``1\,\mathrm{s}`` record length of microseismic data (Figure [#acq1b]) to jointly estimate the location and source time function of both the sources. ### Experiment 1—True velocity experiment In the first experiment, we assume that we have access to the true velocity model. To compare how proposed method works in comparison to Linearized Bregman algoirthm for our problem, we perform 50 iterations of both the algorithms for same value of ``\mu``. Figure [#acq1c] & [#acq1d] are intensity plots corresponding to the linearized Bregman and the proposed method based on L-BFGS. The linearized Bregman based method is able to estimate only one source location correctly, namely at (``1.625\, \mathrm{km}`` , ``0.73\, \mathrm{km}``) (Figure [#acq1c]), whereas the proposed method estimates both source locations correctly after the same number of iterations (Figure [#acq1d]). This clearly demonstrates the faster convergence of the proposed method with compared to the linearized Bregman based method. Also, the proposed method is able to give a good estimation of the source time function at both locations, whereas the linearized Bregman based results only recovers one source (Figure [#Wavelet_t1] & [#Wavelet_t2]). The proposed method is able to preserve the frequency content of the source time function at both the source locations. This is evident from the amplitude spectrum plot of source time function at both the locations (Figure [#Spectrum_t1] & [#Spectrum_t2]). This is because we are working in ``2.5``D instead of ``2``D, which prevents any kind of frequency loss due to ``2``D wave modeling [@song1995freq]. #### Figure: {#acq1} ![](Figure/SEG_Figures/True_Vel_Marm.png){#acq1a width=50%} ![](Figure/SEG_Figures/Shot_true_Marm.png){#acq1b width=50%}\ ![](Figure/SEG_Figures/Est_src_lbr_True_Marm_zoom.png){#acq1c width=50%} ![](Figure/SEG_Figures/Est_src_lbf_True_Marm_zoom.png){#acq1d 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}``. Two black dots indicate the location of two microseismic sources. (b) Microseismic data generated by two microseismic sources activating at different times and with dominant frequencies of ``30.0\, \mathrm{Hz}`` and ``25.0\, \mathrm{Hz}`` respectively. Source location estimated by (c) Linearized Bregman and (d) L-BFGS as the outlier in the zoomed intensity plot. Black cross indicates the true source location #### Figure: {#Wavelet_true} ![](Figure/SEG_Figures/Wavelet_comparison_loc1.png){#Wavelet_t1 width=50%} ![](Figure/SEG_Figures/Wavelet_comparison_loc2.png){#Wavelet_t2 width=50%}\ ![](Figure/SEG_Figures/Spectrum_comparison_loc1.png){#Spectrum_t1 width=50%} ![](Figure/SEG_Figures/Spectrum_comparison_loc2.png){#Spectrum_t2 width=50%} : Comparison of the true wavelet (solid green) with wavelets estimated by Linearized Bregman based method (solid blue) and L-BFGS based method (solid magenta) using true velocity and wavelets estimated by Linearized Bregman based method (dashed blue) and L-BFGS based method (dashed magenta) using smooth velocity at (a) source ``1`` located at (``1.625\, \mathrm{km}`` , ``0.73\, \mathrm{km}``) and (b) source 2 located at (``1.665\, \mathrm{km}`` , ``0.72\, \mathrm{km}``). Comparison of spectrum of true wavelet (solid green) with spectrum of wavelets estimated by Linearized Bregman based method (solid blue) and L-BFGS based method (solid magenta) using true velocity and spectrum of wavelets estimated by Linearized Bregman based method (dashed blue) and L-BFGS based method (dashed magenta) using smooth velocity at (c) source ``1`` located at (``1.625\, \mathrm{km}`` , ``0.73\, \mathrm{km}``) and (d) source 2 located at (``1.665\, \mathrm{km}`` , ``0.72\, \mathrm{km}``) ### Experiment 2—Smooth velocity experiment In the 2nd experiment we assume that we only have access to a smooth velocity model. We still use the data generated in experiment 1 (Figure [#acq1b]), but this time we use a smoothed version of the Marmousi model to jointly estimate the microseismic source location and source time function of both the sources. We use ``125\ \times 125\, \mathrm{m^2}`` ``2``d filter to smooth the Marmousi model. For same value of ``\mu`` , we perform 100 iterations of linearized Bregman based method and the proposed method to jointly estimate location and source-time function of both the microseismic sources. As before, with the linearized Bregman based method we are able to recover only one source location (Figure [#loca]) after 100 iterations whereas with the proposed method, we are able to estimate both the source locations (Figure [#locb]). We get good estimate of the source time function at both the locations with the proposed method, whereas linearized Bregman based method estimates source time function at one source location (Figure [#Wavelet_t1]). The proposed method also preserves the frequency content of the source time function at both the source locations as evident from the amplitude spectrum plot of source time function at both the locations (Figure [#Spectrum_t1] & [#Spectrum_t2]). We observe some noise in the spectrum of estimated source-time function. This is mainly because of the cross-talk between sources within half a wavelength. We can reduce this cross talk by spectral smoothing. With velocity update, we can get further focusing of source field at the true source location, which can further reduce this cross talk. Velocity update is out of scope for this paper but the comparison of results of experiment ``1`` with true velocity and experiment ``2`` with smooth velocity is encouraging. #### Figure: {#loc} ![](Figure/SEG_Figures/Est_src_lbr_Smooth_Marm_zoom.png){#loca width=50%} ![](Figure/SEG_Figures/Est_src_lbf_Smooth_Marm_zoom.png){#locb width=50%} : Source location estimated by (a) Linearized Bregman and (b) L-BFGS as the outlier in zoomed intensity plot.Black cross indicates the true source location \vspace*{-0.5cm} ## Conclusions \vspace*{-0.25cm} We have proposed a method to jointly estimate closely spaced microseismic source locations and their source time function functions. Linearized Bregman based method requires lot of iterations for closely spaced sources. However, L-BFGS based method can locate sources within half a wavelength along with their source time function in lesser iterations. The L-BFGS based method works with dual variable which is cheaper to store. Therefore, we can store dual variable updates without any addition memory storage burden to get better L-BFGS convergence. We have extended our method in ``2.5``D. Lastly, we can estimate location and source time function of more than one sources without any assumptions on number of sources. \vspace*{-0.5cm} ## Acknowledgement \vspace*{-0.25cm} This research was carried out as part of the SINBAD 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}}} ```