Metric operations

Building a metric that is suitable for remeshing may require combining information from different sources or imposing constraints (e.g. minimum or maximum edge lengths, smooth variation with respect to the input mesh, ... ).

Metric Interpolation

In practice, depending on the CFD solver, the metric is either stored on the vertex or the cell center. The definition of an interpolation procedure on metrics is therefore mandatory to be able to compute the metric at any point of the domain. When more than two discrete metrics needs to be interpolated, it is important to use a consistent interpolation framework. For instance, the interpolation needs to be commutative (i.e. the resulting metric does not depend on the order of the interpolation operations between metrics). To that end, the log-Euclidean framework introduced by Arsigny et al. (2006) is used.

Log-Euclidean framework

The authors first define the notion of metric logarithm and matrix exponential. The metric logarithm is defined on the set of metric tensors. For metric tensor $\mathcal{M}=\mathcal{R}\Lambda {\mathcal{R}}^{\top }$, it is given by:

$$\ln (\mathcal{M}):=\mathcal{R}\ln (\Lambda ){\mathcal{R}}^{\top }$$

where $\ln (\Lambda )=\text{diag}(\ln ({\lambda }_i))$. The matrix exponential is defined on the set of symmetric matrices. For any symmetric matrix $S=\mathcal{Q}\Sigma {\mathcal{Q}}^{\top }$, it is given by:

$$\exp (S):=\mathcal{Q}\exp (\Sigma ){\mathcal{Q}}^{\top }$$

where $\exp (\Sigma )=\text{diag}(\exp ({\xi }_i))$. We can now define the logarithmic addition $\oplus$ and the logarithmic scalar multiplication $\odot$:

$${\mathcal{M}}_1\oplus {\mathcal{M}}_2:=\exp (\ln ({\mathcal{M}}_1)+\ln ({\mathcal{M}}_2))$$

$$\alpha \odot \mathcal{M}:=\exp (\alpha \cdot \ln (\mathcal{M}))={\mathcal{M}}^{\alpha }$$

The logarithmic addition is commutative and coincides with matrix multiplication whenever the two tensors ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$ commute in the matrix sense.

Metric interpolation in the log-Euclidean framework

Let $({\mathbf{x}}_i{)}_{i=1,...k}$ be a set of cell centers and $(\mathcal{M}({\mathbf{x}}_i){)}_{i=1,...k}$ their associated metrics. Then, for a point $\mathbf{x}$ of the domain such that:

$$\mathbf{x}=\sum _{i=1}^k{\alpha }_i\cdot {\mathbf{x}}_i\text{ with }\sum _{i=1}^k{\alpha }_i=1$$

the interpolated metric is defined by:

$$\mathcal{M}(\mathbf{x})=\underset{i=1}{\overset{k}{⨁}}{\alpha }_i\odot \mathcal{M}({\mathbf{x}}_i)=\exp \left(\sum _{i=1}^k{\alpha }_i\ln (\mathcal{M}({\mathbf{x}}_i))\right)$$

This interpolation is commutative, but its bottleneck is to perform $k$ diagonalizations and to request the use of the logarithm and the exponential functions which are CPU consuming. However, this procedure is essential to define continuously the metric map on the entire domain. Also, since CODA is a cell centered solver and Tucanos is vertex based, this interpolation definition is crucial to get the metric field at the vertexes of the mesh. Moreover, it has been demonstrated (Arsigny et al., 2006) that this interpolation preserves the maximum principle, i.e., for an edge $\mathbf{ab}$ with endpoint metrics $\mathcal{M}(\mathbf{a})$ and $\mathcal{M}(\mathbf{b})$ such that $\text{det}(\mathcal{M}(\mathbf{a}))<\text{det}(\mathcal{M}(\mathbf{b}))$ then we have:

$$\text{det}(\mathcal{M}(\mathbf{a}))<\text{det}(\mathcal{M}(\mathbf{a}+t\mathbf{ab}))<\text{det}(\mathcal{M}(\mathbf{b})),\forall t\in [0,1]$$

Metric Intersection

In metric-based adaptation, the metric field computed from the error estimation does not take into account geometric or gradation constraints, that would be defined using different metric fields. Metric intersection is used to get a unique metric field. It consists in keeping the most restrictive size constraint in all directions imposed by this set of metrics.

Figure 1: Left, view illustrating the metric intersection procedure with the simultaneous reduction in three dimensions. In red, the resulting metric of the intersection of the blue and green metrics. Right, metric interpolation along a segment where the endpoints metrics are the blue and violet ones. (Loseille 2008)

Formally speaking, let ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$ be two metric tensors given at a point. The metric tensor ${\mathcal{M}}_{1\cap 2}$ corresponding to the intersection of ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$ is the one prescribing the largest possible size under the constraint that the size in each direction is always smaller than the sizes prescribed by ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$. Let us give a geometric interpretation of this operator. Metric tensors are geometrically represented by an ellipse in 2D and an ellipsoid in 3D. But the intersection between two metrics is not directly the intersection between two ellipsoids as their geometric intersection is not an ellipsoid. Therefore, we seek for the largest ellipsoid representing ${\mathcal{M}}_{1\cap 2}$ included in the geometric intersection of the ellipsoids associated with ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$, cf. the Figure above (left). The ellipsoid (metric) verifying this property is obtained by using the simultaneous reduction of two metrics.

Simultaneous reduction

The simultaneous reduction enables to find a common basis $({\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3)$ such that ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$ are congruent to a diagonal matrix in this basis, and then to deduce the intersected metric. To do so, the matrix $\mathcal{N}={\mathcal{M}}_1^{-1}{\mathcal{M}}_2$ is introduced. $\mathcal{N}$ is diagonalizable with real-eigenvalues. The normalized eigenvectors of $\mathcal{N}$ denoted by ${\mathbf{e}}_1$, ${\mathbf{e}}_2$ and ${\mathbf{e}}_3$ constitute a common diagonalization basis for ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$. The entries of the diagonal matrices, that are associated with the metrics ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$ in this basis, are obtained with the Rayleigh formula:

$${\lambda }_i={\mathbf{e}}_i^{\top }{\mathcal{M}}_1{\mathbf{e}}_i\text{ and }{\mu }_i={\mathbf{e}}_i^{\top }{\mathcal{M}}_2{\mathbf{e}}_i$$

Remark: ${\lambda }_i$ and ${\mu }_i$ are not the eigenvalues of ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$. They are the spectral values associated with basis $({\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3)$.

Let $\mathcal{P}=({\mathbf{e}}_1{\mathbf{e}}_2{\mathbf{e}}_3)$ be the matrix with the columns that are the eigenvectors $\{{\mathbf{e}}_i{\}}_{i=\mathrm{1...3}}$ of $\mathcal{N}$. $\mathcal{P}$ is invertible as $({\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3)$ is a basis of ${\mathbb{R}}^3$. We have:

$${\mathcal{M}}_1={\mathcal{P}}^{-\top }\left(\begin{array}{ccc}{\lambda }_1& 0& 0\\ 0& {\lambda }_2& 0\\ 0& 0& {\lambda }_3\end{array}\right){\mathcal{P}}^{-1}\text{ and }{\mathcal{M}}_2={\mathcal{P}}^{-\top }\left(\begin{array}{ccc}{\mu }_1& 0& 0\\ 0& {\mu }_2& 0\\ 0& 0& {\mu }_3\end{array}\right){\mathcal{P}}^{-1}$$

Computing the metric intersection

The resulting intersected metric ${\mathcal{M}}_{1\cap 2}$ is then analytically given by:

$${\mathcal{M}}_{1\cap 2}={\mathcal{M}}_1\cap {\mathcal{M}}_2={\mathcal{P}}^{-\top }\left(\begin{array}{ccc}\max ({\lambda }_1,{\mu }_1)& 0& 0\\ 0& \max ({\lambda }_2,{\mu }_2)& 0\\ 0& 0& \max ({\lambda }_3,{\mu }_3)\end{array}\right){\mathcal{P}}^{-1}$$

The ellipsoid associated with ${\mathcal{M}}_{1\cap 2}$ is the largest ellipsoid included in the geometric intersection region of the ellipsoids associated with ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$.

Numerically, to compute ${\mathcal{M}}_{1\cap 2}$, the real-eigenvalues of $\mathcal{N}$ are first evaluated with a Newton algorithm. Then, the eigenvectors of $\mathcal{N}$, which define $\mathcal{P}$, are computed using the algebra notions of image and kernel spaces.

See proof

Metric gradation

Metric fields may have huge variations or may be quite irregular when evaluated from numerical solutions that present discontinuities or steep gradients. This makes the generation of a unit mesh difficult or impossible, thus leading to poor quality anisotropic meshes. Generating high-quality anisotropic meshes requires to smooth the metric field by bounding its variations in all directions. It also helps with the convergence of the implicit solver. Mesh gradation (Alauzet, 2010) consists in reducing in all directions the size prescribed at any points if the variation of the metric field is larger than a fixed threshold.

Spanning a metric field

Let $\mathbf{p}$ be a point of a domain $\Omega$ supplied with a metric ${\mathcal{M}}_{\mathbf{p}}$ and $\beta$ the specified gradation. Two laws governing the metric growth in the domain are proposed. In the first one, the metric growth is homogeneous in the Euclidean metric field defined by ${\mathcal{M}}_{\mathbf{p}}$. The second metric growth is homogeneous in the physical space, i.e., the classical Euclidean space.

The first law associates for any point $\mathbf{x}$ of the domain a unique scale factor with the metric given by:

$${\eta }^2(\mathbf{px})=(1+l_p(\mathbf{px})\cdot \ln (\beta ){)}^{-2}=(1+\sqrt{{\mathbf{px}}^{\top }{\mathcal{M}}_{\mathbf{p}}\mathbf{px}}\cdot \ln (\beta ){)}^{-2}$$

With this formulation, each pointwise metric ${\mathcal{M}}_{\mathbf{p}}$ spans a global continuous smooth metric field all over domain $\Omega$ parametrized by the given gradation value $\beta$:

$$({\mathcal{M}}_{\mathbf{p}}(\mathbf{x}){)}_{\mathbf{x}\in \Omega }\text{with }{\mathcal{M}}_{\mathbf{p}}(\mathbf{x})={\eta }^2(\mathbf{px}){\mathcal{M}}_{\mathbf{p}}$$

In this case, the resulting metric field grows homogeneously in the Euclidean metric space defined by ${\mathcal{M}}_{\mathbf{p}}$ as the scale factor depends on the length of segment $\mathbf{px}$ with respect to ${\mathcal{M}}_{\mathbf{p}}$. As a result, the shape of the metric is kept unchanged while growing. This law conserves the same anisotropic ratio.

For the second law, we associate independently a growth factor with each eigenvalue of ${\mathcal{M}}_{\mathbf{p}}=\mathcal{R}\Lambda {\mathcal{R}}^{\top }$, with $\Lambda =\text{diag}({\lambda }_i{)}_{i=1,3}$:

$${\eta }_i^2(\mathbf{px})=(1+\sqrt{{\lambda }_i}||\mathbf{px}|{|}_2\cdot \ln (\beta ){)}^{-2}$$

the grown metric at $\mathbf{x}$ is given by:

$${\mathcal{M}}_{\mathbf{p}}(\mathbf{x})=\mathcal{RN}(\mathbf{px})\Lambda {\mathcal{R}}^{\top }\text{where }\mathcal{N}(\mathbf{px})=\left(\begin{array}{ccc}{\eta }_1^2(\mathbf{px})& 0& 0\\ 0& {\eta }_2^2(\mathbf{px})& 0\\ 0& 0& {\eta }_3^2(\mathbf{px})\end{array}\right)$$

The resulting metric field grows homogeneously in the physical space. Indeed, each eigenvalue grows similarly in all directions, as the factor ${\eta }_i$ depends only on the distance (in the physical space) from the original point. Consequently, the shape, i.e., the anisotropic ratio of the metric, is no more preserved as the eigenvalues are growing separately and differently. This law gradually makes the metric more and more isotropic as it gradually propagates in the domain.

Alauzet et al. (2021) suggest to mix these two laws to achieve an efficient metric gradation algorithm. For this new law, a growth factor is associated independently with each eigenvalue of ${\mathcal{M}}_{\mathbf{p}}$:

$${\eta }_i^2(\mathbf{px})=((1+\sqrt{{\lambda }_i}||\mathbf{px}|{|}_2\cdot \ln (\beta ){)}^t(1+l_p(\mathbf{px})\cdot \ln (\beta ){)}^{1-t}{)}^{-2}$$

The authors consider $t=\frac{1}{8}$ within their numerical examples. For the applications in the work, we will use $t=1.$

Metric reduction

The reduced metric at a point $\mathbf{x}$ of the domain $\Omega$ is given by the strongest size constraint imposed by the metric at $\mathbf{x}$ and by the spanned metrics (parametrized by the given size gradation) of all the other points of the domain at $\mathbf{x}$:

$$\widetilde{\mathcal{M}(\mathbf{x})}=\left(\bigcap _{\mathbf{p}\in \Omega }{\mathcal{M}}_{\mathbf{p}}(\mathbf{x})\right)\cap \mathcal{M}(\mathbf{x})$$

Practical implementation are given in Alauzet (2010).

Thresholds on the local sizes

In order to ensure that the metric field will generate meshes with edge sizes between $h_{min}$ and $h_{max}$, the metric $\mathcal{M}=R^{T}diag(s_1,\cdots s_d)R$ with $R=({\mathbf{v}}_1|\cdots |{\mathbf{v}}_d)$ is modified as $$\mathcal{T}(\mathcal{M},h_{min},h_{max})=R^{T}diag(s_1,\cdots s_d)R$$ with ${\tilde{s}}_i=\min \left(\max \left(s_i,\frac{1}{h_{max}^2}\right),\frac{1}{h_{min}^2}\right)$.

Element-implied metric

Given an element $K$, defined by its vertices $({\mathbf{x}}_0,\cdots ,{\mathbf{x}}_d)$, the element-implied metric is the metric for which the element is equilateral

If $J_{K,\Delta }$ is the jacobian of the transformation from the reference equilateral element to the physical element, the element-implied metric is $${\mathcal{M}}_K=(J_{K,\Delta }{J}_{K,\Delta }^T{)}^{-1}$$

In practice $J_{K,\Delta }$ is computed as $J_{K,\Delta }=J_{K,\perp }{J}_{\perp ,\Delta }$ where - $J_{K,\perp }=({\mathbf{x}}_1-{\mathbf{x}}_0|\cdots |{\mathbf{x}}_d-{\mathbf{x}}_0)$ is the jacobian of transformation from the reference orthogonal element to the physical element - $J_{\perp ,\Delta }$ is the jacobian of transformation from the reference equilateral element to the reference orthogonal element (independent of $K$ )

Controlling the step between two metrics

In the remeshing process, it may be interesting to control the step between the element-implied and target metrics. Such control is available within the remesher Tucanos and specified by the parameter $f$.

Given two metric fields ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$, the objective is to find $\mathcal{M}=\mathcal{L}({\mathcal{M}}_1,{\mathcal{M}}_2,f)$ as close as possible from ${\mathcal{M}}_2$ such that, for all edges $e$:

$$\frac{1}{f}\le \sqrt{\frac{e^T\mathcal{M}e}{e^T{\mathcal{M}}_{1}e}}\le f$$

i.e. to have:

$$\frac{1}{f}\le {\lambda }_{\min }({\mathcal{M}}_1^{-\frac{1}{2}}\mathcal{M}{\mathcal{M}}_1^{-\frac{1}{2}})\le {\lambda }_{\max }({\mathcal{M}}_1^{-\frac{1}{2}}\mathcal{M}{\mathcal{M}}_1^{-\frac{1}{2}})\le f$$

Practically, we define as close as possible as minimizing the Froebenius norm $\Vert {\mathcal{M}}_1^{-\frac{1}{2}}(\mathcal{M}-{\mathcal{M}}_2){\mathcal{M}}_1^{-\frac{1}{2}}{\Vert }_F$. The optimal $\mathcal{M}$ is then computed as follows:

• compute $\mathcal{N}:={\mathcal{M}}_1^{-\frac{1}{2}}{\mathcal{M}}_2{\mathcal{M}}_1^{-\frac{1}{2}}$,
• compute the eigenvalue decomposition $QD{Q}^T=\mathcal{N}$, with $D=\text{diag}({\lambda }_i)$
• compute ${\mathcal{N}}^{\ast }:=Q\text{diag}({\hat{\lambda }}_i)Q^T$ where ${\hat{\lambda }}_i:=\min (\max ({\lambda }_i,\frac{1}{f}),f)$,
• compute ${\mathcal{M}}^{\ast }:={\mathcal{M}}_1^{\frac{1}{2}}{\mathcal{N}}^{\ast }{\mathcal{M}}_1^{\frac{1}{2}}$.

Figure 2: Representation of the controlled metric $\mathcal{L}({\mathcal{M}}_1,{\mathcal{M}}_2,f)$ by stepping with stepping parameter $f=2$. The thin blue lines represent $2{\mathcal{M}}_1$ and $\frac{1}{2}{\mathcal{M}}_1$, X. Garnaud

See proof

Metric scaling

• The "ideal" metric $\mathcal{M}$ is modified using a series of constraints
  • Minimum and maximum edge sizes $h_{min}$ and $h_{max}$,
  • Maximum gradation $\beta$,
  • Maximum step $f$ with respect to the element implied metric ${\mathcal{M}}_i$,
  • A minimum curvature radius to edge lenth ratio $r$ (and optionally $y_{target}^{+}$), giving a fixed metric ${\mathcal{M}}_f$
• A scaling $\alpha$ is applied, such that the actual metric used for remeshing is $$\tilde{\mathcal{M}}=\mathcal{L}(\mathcal{T}(\alpha \mathcal{M},h_{min},h_{max})\cap {\mathcal{M}}_f,{\mathcal{M}}_i,f)$$
• $\alpha$ is chosen to have a given complexity $$\mathcal{C}(\tilde{\mathcal{M}})=\mathcal{N}$$
• A simple bisection is used to compute $\alpha$ if all the constraints can be met.

Maths

Notations

• One denotes the Frobenius norm of a matrix $A\in {\mathbb{R}}^{n\times n}$ by $\Vert A{\Vert }_F:=\mathrm{tr}(A^{T}A{)}^{1/2}=({\sum }_{i,j}{A}_{i,j}^2{)}^{1/2}$. This norm is associated with the scalar product $(A,B{)}_F:=\mathrm{tr}(A^{T}B)$.
• One denotes by $\le$ the partial order in the symmetric matrices space ${\mathcal{S}}_n(\mathbb{R})$: $S_1\le S_2$ if $S_2-S_1\in {\mathcal{S}}_n^{+}(\mathbb{R})$, i.e. for all $x\in {\mathbb{R}}^n$, $x^T{S}_{1}x\le x^T{S}_{2}x$.

Projection on a closed convex subset of the symmetric matrices

The following theorem has been proven in Computing a nearest symmetric positive semidefinite matrix, Higham Nicholas J, 1988:

Theorem 1 Let ${\mathcal{S}}_{\delta }:=\{S\in {\mathcal{S}}_n({\mathbb{R}}^n):\delta I_n\le S\}$. Then, ${\mathcal{S}}_{\delta }$ is a closed convex subset of the space ${\mathcal{S}}_n({\mathbb{R}}^n)$ and then there exists a projection application $P_{\delta }:{\mathcal{S}}_n({\mathbb{R}}^n)\to {\mathcal{S}}_{\delta }$ such that $$\forall A\in {\mathcal{S}}_n(\mathbb{R}),\Vert A-P_{\delta }(A){\Vert }_F=\inf \{\Vert A-S{\Vert }_F:S\in {\mathcal{S}}_{\delta }\}.>$$ Let $A\in {\mathcal{S}}_n(\mathbb{R})$ and $A=Q^{T}DQ$ its diagonalisation in an orthonormal basis: $Q\in {\mathcal{O}}_n(\mathbb{R})$, $D=\mathrm{diag}({\lambda }_i)$. Then, $$P_{\delta }(A)=Q\hat{D}{Q}^T,\text{where}\hat{D}:=\mathrm{diag}(\max ({\lambda }_i,\delta )).$$

Proof

Let $A\in {\mathcal{S}}_n(\mathbb{R})$ and its diagonalisation in a orthonomal basis $P^{T}AP=D=\mathrm{diag}({\lambda }_i)$. As $\Vert .{\Vert }_F$ is invariant by linear isometry ($\Vert PAQ{\Vert }_F^2=\mathrm{tr}((PAQ{)}^{T}PAQ)=\mathrm{tr}(Q^T{A}^T{P}^{T}PAQ)=\mathrm{tr}(Q{Q}^T{A}^{T}A)=\Vert A{\Vert }_F^2$, one gets $\Vert A-S{\Vert }_F^2=\Vert D-Q^{T}SQ{\Vert }_F^2$ and then the change of variable $S\in {\mathcal{S}}_n(\mathbb{R})\mapsto X=Q^{T}SQ\in {\mathcal{S}}_n(\mathbb{R})$ together with $\delta I\le S\iff \delta I\le Q^{T}SQ$, leads to the following problem $$\inf \{\Vert D-X{\Vert }_F^2:X\in {\mathcal{S}}_n(\mathbb{R}),\delta I\le X\}.$$ This problem can be computed explicitly. Since $X_{ii}=e_i^{T}X{e}_i\ge \delta e_i^T{e}_i=\delta$, ($e_i=({\delta }_{ij}{)}_j$ the canonical basis vectors) $$\Vert D-X{\Vert }_F^2=\sum _{i=1}^n({\lambda }_i-X_{ii}{)}^2+2\sum _{i<j}{X}_{ij}^2\ge \sum _{i=1}^n({\lambda }_i-X_{ii}{)}^2\ge \sum _{i:{\lambda }_i<\delta }({\lambda }_i-\delta {)}^2.$$ This lower bound is reached for the following diagonal matrix $X:=\hat{D}:=\mathrm{diag}({\hat{\lambda }}_i)$, where ${\hat{\lambda }}_i:=\max ({\lambda }_i,\delta )$, and moreover, $\hat{D}\ge \delta I$ and is in ${\mathcal{S}}_n(\mathbb{R})$, which ends the proof.

Generalization

One can generalize the Theorem above as follow:

Theorem 2 Given two symmetric real numbers $0<\alpha <\beta$, let ${\mathcal{S}}_{\alpha ,\beta }:=\{S\in {\mathcal{S}}_n(\mathbb{R}):\alpha I_n\le S\le \beta I_n\}$. ${\mathcal{S}}_{\alpha ,\beta }$ is a compact convex subset of ${\mathcal{S}}_n(\mathbb{R})$ and then let us denote $P_{\alpha ,\beta }:{\mathcal{S}}_n(\mathbb{R})\to {\mathcal{S}}_{\alpha ,\beta }$ the associated projection application: $$\forall A\in {\mathcal{S}}_n(\mathbb{R}),\Vert A-P_{\alpha ,\beta }(A){\Vert }_F=\inf \{\Vert A-S{\Vert }_F:S\in {\mathcal{S}}_{\alpha ,\beta }\}.>$$ Given $A\in {\mathcal{S}}_n(\mathbb{R})$ and $A=QD{Q}^T$ its diagonalisation in an orthonormal basis: $Q\in {\mathcal{O}}_n(\mathbb{R})$, $D=\mathrm{diag}({\lambda }_i)$. Then, $$P_{\alpha ,\beta }(A)=Q\hat{D}{Q}^T,\text{where}\hat{D}:=\mathrm{diag}({\hat{\lambda }}_i),{\hat{\lambda }}_i:=\min \left(\max \right({\lambda }_i,\alpha ),\beta ).$$

Proof

The proof follows the one above very closely. The lower bound is got using the inequality: $\alpha e_i^T{I}_n{e}_i\le e_i^{T}X{e}_i\le \beta e_i^T{I}_n{e}_i\Rightarrow \alpha \le X_{ii}\le \beta$: $$\Vert D-X{\Vert }_F^2\ge \sum _{i=1}^n({\lambda }_i-X_{ii}{)}^2+2\sum _{i<j}{X}_{ij}^2$$ $$\ge \sum _{i=1}^n({\lambda }_i-X_{ii}{)}^2$$ $$\ge \sum _{i:{\lambda }_i<\alpha }({\lambda }_i-\alpha {)}^2+\sum _{i:{\lambda }_i>\beta }({\lambda }_i-\beta {)}^2,$$ which is reached for $\hat{D}:=\mathrm{diag}({\hat{\lambda }}_i)$ where ${\hat{\lambda }}_i=\alpha$ if ${\lambda }_i<\alpha$, ${\hat{\lambda }}_i=\beta$ if ${\lambda }_i>\beta$ and otherwise ${\hat{\lambda }}_i={\lambda }_i$. Again, the matrix $\hat{D}$ is symmetric and indeed in ${\mathcal{S}}_{\alpha ,\beta }$ (because since $\hat{D}$ is diagonal, $\alpha \le {\hat{\lambda }}_i\le \beta \iff \alpha I_n\le \hat{D}\le \beta I_n$) and thus is the unique solution of the problem (existence and uniqueness being given by the projection on closed convex sets in a Hilbert space).

Metric optimization problem

Given two metrics ${\mathcal{M}}_1,{\mathcal{M}}_2\in {\mathcal{S}}_n^{++}(\mathbb{R})$. The deformation ratio between a metric $\mathcal{M}\in S_n^{++}(\mathbb{R})$ and ${\mathcal{M}}_0$ is defined as a function of x∈Rn</span>{0} r(M,M0​;x):=(xTM0​xxTMx​)1/2.

The problem.

Given two metrics ${\mathcal{M}}_1,{\mathcal{M}}_2\in {\mathcal{S}}_n^{++}(\mathbb{R})$, and two real numbers $0<\alpha <\beta$, let ${\mathcal{R}}_{\alpha ,\beta }({\mathcal{M}}_0)$ be the set of metrics which are not to stretched with regard to ${\mathcal{M}}_0$: Rα,β​(M0​):={M∈Sn​(R):∀x∈Rn</span>{0},α≤r(M,M0​;x)2≤β}. Again, this set is compact and convex. Given a matrix norm $$|||.|||$$, we want to find $${\mathcal{M}}^{\ast }\in {\mathcal{R}}_{\alpha ,\beta }({\mathcal{M}}_0)$$ which minimize the distance to $${\mathcal{M}}_1$$, i.e. find $${\mathcal{M}}^{\ast }\in {\mathcal{R}}_{\alpha ,\beta }({\mathcal{M}}_0)$$ such that: $$|||{\mathcal{M}}_1-{\mathcal{M}}^{\ast }|||=\inf \{|||{\mathcal{M}}_1-\mathcal{M}|||:\mathcal{M}\in {\mathcal{R}}_{\alpha ,\beta }({\mathcal{M}}_0)\}$$

Lemma : Let $\mathcal{M},{\mathcal{M}}_0\in {\mathcal{S}}_n^{++}(\mathbb{R})$ two metrics > and $\alpha ,\beta \in \mathbb{R}$. The function x∈Rn</span>{0}↦r(M,M0​;x)2 is bounded and reaches its bounds which are the minimum and maximum of the eigenvalues of M0>​−1/2MM0−1/2​. Thus, >∀x∈Rn</span>{0},>α≤r(M,M0​;x)≤β> ⇔>αIn​≤M0−1/2​MM0−1/2​≤β>In​> ⇔>αM0​≤M≤βM0​>

Thanks to this lemma, the admissible set ${\mathcal{R}}_{\alpha ,\beta }({\mathcal{M}}_0)$ can be re-written as $${\mathcal{R}}_{\alpha ,\beta }({\mathcal{M}}_0)=\{\mathcal{M}\in {\mathcal{S}}_n(\mathbb{R}):\alpha {\mathcal{M}}_0\le \mathcal{M}\le \beta {\mathcal{M}}_0\}$$ $$=\{\mathcal{M}\in {\mathcal{S}}_n(\mathbb{R}):\alpha I_n\le {\mathcal{M}}_0^{-1/2}\mathcal{M}{\mathcal{M}}_0^{-1/2}\le \beta I_n\}.$$

In the ${\mathcal{M}}_0$-twisted Frobenius norm

First, one studies the case of choosing $$|||\mathcal{M}|||:=\Vert {\mathcal{M}}_1-\mathcal{M}{\Vert }_{{\mathcal{M}}_0,F}:=\Vert {\mathcal{M}}_0^{-1/2}\mathcal{M}{\mathcal{M}}_0^{-1/2}{\Vert }_F.$$ The above problem reads as follows. Find ${\mathcal{M}}^{\ast }\in {\mathcal{R}}_{\alpha ,\beta }({\mathcal{M}}_0)$ such that: $$|||{\mathcal{M}}_1-{\mathcal{M}}^{\ast }|||=\inf \{||{\mathcal{M}}_0^{-1/2}({\mathcal{M}}_1-\mathcal{M}){\mathcal{M}}_0^{-1/2}|{|}_F:\mathcal{M}\in {\mathcal{S}}_n(\mathbb{R}),{\mathcal{M}}_0^{-1/2}\mathcal{M}{\mathcal{M}}_0^{-1/2}\in {\mathcal{S}}_{\alpha ,\beta }\}.$$ It is then clear that the change of variable $\mathcal{M}\in {\mathcal{S}}_n(\mathbb{R})\to \mathcal{N}:={\mathcal{M}}_0^{-1/2}\mathcal{M}{\mathcal{M}}_0^{-1/2}\in {\mathcal{S}}_n(\mathbb{R})$ leads to a problem of the form of Theorem 2. Denoting ${\mathcal{N}}_1:={\mathcal{M}}_0^{-1/2}{\mathcal{M}}_1{\mathcal{M}}_0^{-1/2}$ the problem reads find ${\mathcal{M}}^{\ast }={\mathcal{M}}_0^{1/2}{\mathcal{N}}^{\ast }{\mathcal{M}}_0^{1/2}\in {\mathcal{R}}_{\alpha ,\beta }({\mathcal{M}}_0)$ such that $$\Vert {\mathcal{N}}_1-{\mathcal{N}}^{\ast }{\Vert }_F=\inf \{\Vert {\mathcal{N}}_1-\mathcal{N}{\Vert }_F:\mathcal{N}\in {\mathcal{S}}_{\alpha ,\beta }\}.$$ Then, from Theorem 2, the solution is $${\mathcal{M}}^{\ast }={\mathcal{M}}_0^{1/2}{P}_{\alpha ,\beta }({\mathcal{N}}_1){\mathcal{M}}_0^{1/2}$$ and can be computed as follows:

• Compute ${\mathcal{N}}_1:={\mathcal{M}}_0^{-1/2}{\mathcal{M}}_1{\mathcal{M}}_0^{-1/2}$,
• Find $Q\in {\mathcal{O}}_n(\mathbb{R})\text{ and }D=\mathrm{diag}({\lambda }_i)$ such that $QD{Q}^T={\mathcal{N}}_1$.
• Compute ${\mathcal{N}}^{\ast }:=P_{\alpha ,\beta }({\mathcal{N}}_1)=Q\mathrm{diag}({\hat{\lambda }}_i)Q^T,\text{where }{\hat{\lambda }}_i:=\min \left(\max \right({\lambda }_i,\alpha ),\beta )$.
• Compute ${\mathcal{M}}^{\ast }:={\mathcal{M}}_0^{1/2}{\mathcal{N}}^{\ast }{\mathcal{M}}_0^{1/2}$

Remark: Maybe a more suitable formulation would be minimizing the Frobenius norm $|||{\mathcal{M}}_1-\mathcal{M}|||:=\Vert {\mathcal{M}}_1-\mathcal{M}{\Vert }_F$. This leads to a different problem. For instance, if $${\mathcal{M}}_0=\left(\begin{array}{cc}1& 0\\ 0& 1/\epsilon \end{array}\right)$$ where $0<\epsilon \ll 1$. One gets, denoting the error $\mathcal{E}:=(e_{ij}):={\mathcal{M}}_1-{\mathcal{M}}^{\ast }$ $$\Vert \mathcal{E}{\Vert }_{{\mathcal{M}}_0,F}^2=e_{11}^2+2\epsilon e_{12}^2+{\epsilon }^2{e}_{22}^2,$$ instead of $$\Vert \mathcal{E}{\Vert }_F^2=e_{11}^2+2e_{12}^2+e_{22}^2.$$ The two distance measures give very different results... does it matter in practice?

Intersection of metrics

Problem.

Given two metrics ${\mathcal{M}}_1,{\mathcal{M}}_2\in {\mathcal{S}}_n^{++}(\mathbb{R})$, we want to find a metric ${\mathcal{M}}_{1\cap 2}\in {\mathcal{S}}_n^{++}(\mathbb{R})$ which respects the most restrictive length of ${\mathcal{M}}_1$ and ${\mathcal{M}}_1$. Moreover, we want that this metric being maximal. More formally, the lenght restriction can be expressed by ∀x∈Rn</span>{0},xTM1∩2​x≤xTM1​xandxTM1∩2​x≤xTM2​x, indeed, $$\sqrt{x^T\mathcal{M}x}$$ is the length of $$x$$ associated with a metric $$\mathcal{M}$$. The local volume associated with a metric $$\mathcal{M}$$ is its determinant: $$\det \mathcal{M}$$, so that the optimization problem reads $$\max \det {\mathcal{M}}_{1\cap 2}$$ such that $${\mathcal{M}}_{1\cap 2}\in {\mathcal{S}}_n(\mathbb{R}),$$ $$0\le {\mathcal{M}}_{1\cap 2}\le {\mathcal{M}}_1,$$ $$0\le {\mathcal{M}}_{1\cap 2}\le {\mathcal{M}}_2.$$

Solution of the problem.

First, one multiplies by ${\mathcal{M}}_1^{-1/2}\in {\mathcal{S}}_n^{++}(\mathbb{R})$ on each side of the constraints, which leads to $$0\le {\mathcal{M}}_1^{-1/2}{\mathcal{M}}_{1\cap 2}{\mathcal{M}}_1^{-1/2}\le I_n$$ $$0\le {\mathcal{M}}_1^{-1/2}{\mathcal{M}}_{1\cap 2}{\mathcal{M}}_1^{-1/2}\le {\mathcal{M}}_1^{-1/2}{\mathcal{M}}_2{\mathcal{M}}_1^{-1/2}.$$ Then, denoting ${\overline{\mathcal{M}}}_2:={\mathcal{M}}_1^{-1/2}{\mathcal{M}}_2{\mathcal{M}}_1^{-1/2}$ and applying the change of variable $\mathcal{M}\mapsto \mathcal{N}:={\mathcal{M}}_1^{-1/2}\mathcal{M}{\mathcal{M}}_1^{-1/2}$, one gets the simplified problem: $$\max \det \left({\mathcal{M}}_1^{1/2}\mathcal{N}{\mathcal{M}}_1^{1/2}\right)$$ such that $$\mathcal{N}\in {\mathcal{S}}_n(\mathbb{R}),$$ $$0\le \mathcal{N}\le I_n,$$ $$0\le \mathcal{N}\le {\overline{\mathcal{M}}}_2.$$

This problem can be even more simplified using $\det \left({\mathcal{M}}_1^{1/2}\mathcal{N}{\mathcal{M}}_1^{1/2}\right)=\det {\mathcal{M}}_1\det \mathcal{N}$ (does not change the $\max$ objective as $\det {\mathcal{M}}_1>0$) and by diagonalization of ${\overline{\mathcal{M}}}_2$ in an orthonomal basis: $${\overline{\mathcal{M}}}_2=P^{T}DP,\text{with }{P}^{T}P=I_n.$$ Indeed, $0\le \mathcal{N}\le P^{T}DP\Rightarrow 0\le P\mathcal{N}{P}^T\le D$ and $0\le \mathcal{N}\le I_n\Rightarrow 0\le P\mathcal{N}{P}^T\le I_n$, and let us change again the variable $\mathcal{X}:=P\mathcal{N}{P}^T$, with $\det \mathcal{N}=\det \mathcal{X}$, the problem is finally: $$\max \det \mathcal{X}$$ such that $$\mathcal{X}\in {\mathcal{S}}_n(\mathbb{R}),$$ $$0\le \mathcal{X}\le I_n,$$ $$0\le \mathcal{X}\le D.$$

This problem can now be solved explicitly by finding an upper bound. Let $\mathcal{Y}\in {\mathcal{S}}_n(\mathbb{R})$ such that $\mathcal{X}={\mathcal{Y}}^2$, then Hadamard inequality provides $\det \mathcal{X}=(\det \mathcal{Y}{)}^2\le {\prod }_j|{\mathcal{Y}}_{\bullet i}{|}_2^2$. Moreover, ${\mathcal{X}}_{jj}={\sum }_k{\mathcal{Y}}_{jk}{\mathcal{Y}}_{kj}=|{\mathcal{Y}}_{\bullet j}{|}_2^2$, from which $\det \mathcal{X}\le {\prod }_i{\mathcal{X}}_{ii}$ and then, using the two constraints: $0\le \mathcal{X}\le I_n\text{ and }0\le \mathcal{X}\le D\Rightarrow 0\le {\mathcal{X}}_{ii}\le \min (D_{ii},1)$, one gets the upper bound $$\det \mathcal{X}\le \prod _{i=1}^n{\mathcal{X}}_{ii}\le \prod _{i=1}^n\min (D_{ii},1),$$ reached for the diagonal matrix $D^{\ast }:=\mathrm{diag}(\min (D_{ii},1))$. Moreover, this matrix is symmetric and verifies $0\le D^{\ast }\le I_n$ and $0\le D^{\ast }\le D$.

Finally, applying the two changes of variables, the unique solution of the problem is $${\mathcal{M}}_{1\cap 2}={\mathcal{M}}_1^{1/2}{P}^T{D}^{\ast }P{\mathcal{M}}_1^{1/2}$$

Hadamard inequality

Given a matrix $M=[v_1|\dots |v_n]\in {\mathbb{C}}^{n\times n}$, one has $$\left|\det M\right|\le \prod _{i=1}^n|v_i{|}_2,$$ with equality if the $(v_1,\dots ,v_n)$ are orthogonal.

The result if clear if $M$ is singular, if not let $\hat{M}:=[{\hat{v}}_1|\dots |{\hat{v}}_n]$, ${\hat{v}}_i=v_i/|v_i{|}_2$ so that one has to show that $|\det \hat{M}|\le 1$. ${\hat{M}}^{\ast }\hat{M}\in {\mathcal{H}}_n(\mathbb{C})$ is hermitian hence diagonalizable in an orthonomal basis, denoting $({\lambda }_1,\dots ,{\lambda }_n)$ its eigenvalues, one has $$|\det \hat{M}{|}^2=\det {\hat{M}}^{\ast }\hat{M}=\prod _{i=1}^n{\lambda }_i$$ $$\le {\left(\frac{1}{n}\sum _{i=1}^n{\lambda }_i\right)}^n={\left(\frac{1}{n}\mathrm{tr}\left({\hat{M}}^{\ast }M\right)\right)}^n={\left(\frac{1}{n}\sum _{i=1}^n\sum _{k=1}^n\overline{{\hat{M}}_{ki}}{\hat{M}}_{ki}\right)}^n$$ $$={\left(\frac{1}{n}\sum _{i=1}^n\sum _{k=1}^n|{\hat{M}}_{ki}{|}^2\right)}^n={\left(\frac{1}{n}\sum _{i=1}^n|{\hat{v}}_i{|}_2^2\right)}^n$$