Simple algorithm for correction of a finite matrix group known approximately

TL;DR

This paper presents a simple iterative correction algorithm for reconstructing finite unitary matrix groups from approximate elements, ensuring consistency with group properties and additional conditions, implemented in open source Python.

Contribution

Introduces a novel correction formula and iterative method for reconstructing finite matrix groups from approximate data, with practical implementation.

Findings

Fast convergence of the correction algorithm.

Effective reconstruction of group elements from approximate data.

Open source Python implementation available.

Abstract

In this paper we have considered a finite unitary matrix group with exact elements being unknown and only approximate elements available. Such a group becomes inconsistent with its own multiplication table. We found simple correction formula for such group. When applied iteratively this formula gives fast convergence of the group elements and allows to perform the group reconstruction. Next we considered small unitary rotation of entire group which makes the group consistent with a set of additional conditions. For example, we might demand that group elements act on certain set of vectors in a predefined manner. Again, iterative procedure based on this correction shows fast convergence. All algorithms developed in this paper were implemented in a Python library which is available as open source software.

Equations

$$\begin{array}{r}\sum _j{F}_{ij}{G}_j^{-1}=N\delta G_i+G_i\sum _j\delta G_j{G}_j^{-1}\\ -\sum _j\delta G_{(ij)}{G}_j^{-1}{G}_i^{-1}{G}_i=N\delta G_i+[G_i,D_1]\end{array}$$

$$\begin{array}{r}\sum _j{G}_j^{-1}{F}_{ji}=N\delta G_i+\sum _j\delta G_j{G}_j^{-1}{G}_i\\ -\sum _j{G}_i{G}_i^{-1}{G}_j^{-1}\delta G_{(ji)}=N\delta G_i+[D_2,G_i]\end{array}$$

$$2N\delta G_i=\sum _j(F_{ij}{G}_j^{-1}+G_j^{-1}{F}_{ji})+[G_i,\sum _j[G_j^{-1},\delta G_j]]$$

$$\delta G_i=(1+R{)}^{+}{G}_i(1+R)-G_i=[G_i,R]$$

$$\delta G_i=\frac{1}{2N}\sum _j\left(F_{ij}(G_j^{\prime }{)}^{-1}+(G_j^{\prime }{)}^{-1}{F}_{ji}\right)$$

$$S_M=\sqrt{\sum _{ij}||F_{ij}|{|}^2}$$

$$S=\sum _{i}Tr\left[(e^{-R}{G}_i{e}^R-Q_i{)}^{+}(e^{-R}{G}_i{e}^R-Q_i)\right]-Tr\left[\lambda (R+R^{+})\right]$$

$$\begin{array}{r}\frac{\partial S}{\partial R}=\sum _i\{2(R+R^{+})-G_i^{+}(R+R^{+})G_i-G_i(R+R^{+})G_i^{+}+[G_i,Q_i^{+}]\\ -\frac{1}{2}\left(R(G_i{Q}_i^{+}+Q_i^{+}{G}_i)+(G_i{Q}_i^{+}+Q_i^{+}{G}_i)R-2G_{i}R{Q}_i^{+}-2Q_i^{+}R{G}_i\right)\}-\lambda =0\end{array}$$

$$\begin{array}{r}\frac{\partial S}{\partial R^{+}}=\sum _i\{2(R+R^{+})-G_i^{+}(R+R^{+})G_i-G_i(R+R^{+})G_i^{+}-[G_i^{+},Q_i]\\ -\frac{1}{2}\left(R^{+}(G_i^{+}{Q}_i+Q_i{G}_i^{+})+(G_i^{+}{Q}_i+Q_i{G}_i^{+})R^{+}-2G_i^{+}{R}^{+}{Q}_i-2Q_i{R}^{+}{G}_i^{+}\right)\}-\lambda =0\end{array}$$

$$\sum _i\left\{R{H}_i+H_{i}R-(G_{i}R{Q}_i^{+}+Q_i^{+}R{G}_i+G_i^{+}R{Q}_i+Q_{i}R{G}_i^{+})\right\}$$

$$=\sum _i[G_i,Q_i^{+}]+[G_i^{+},Q_i]$$

$$H_i=\frac{1}{2}\left(G_i{Q}_i^{+}+Q_i^{+}{G}_i+Q_i{G}_i^{+}+G_i^{+}{Q}_i\right)$$

$$A\otimes B=\left(\begin{array}{ccc}{A}_{11}B& \dots & A_{1n}B\\ 3& & \\ A_{n1}B& \dots & A_{nn}B\end{array}\right)$$

$$P=\sum _i[G_i,Q_i^{+}]+[G_i^{+},Q_i]$$

$$L=\sum _i\left\{H_i^T\otimes I_n+I_n\otimes H_i-(Q_i^{\ast }\otimes G_i+G_i^T\otimes Q_i^{+}+Q_i^T\otimes G_i^{+}+G_i^{\ast }\otimes Q_i)\right\}$$

$$4NR-4\sum _i{G}_{i}R{G}_i^{+}=P$$

$$4NR-4\sum _i{G}_{i}R{G}_i^{+}=\alpha R=(4N-4\mu )R$$

$$R=\frac{1}{4N}\sum _i\left([G_i,Q_i^{+}]+[G_i^{+},Q_i]\right)$$

$$S=\sum _{ij}<G_i{a}_{ij}-b_{ij}|G_i{a}_{ij}-b_{ij}>-Tr\left[\lambda (R+R^{+})\right]$$

$$=\sum _{i}Tr\{G_i\left(\sum _j|a_{ij}><a_{ij}|\right)G_i^{+}-G_i\left(\sum _j|a_{ij}><b_{ij}|\right)$$

$$-\left(\sum _j|b_{ij}><a_{ij}|\right)G_i^{+}+\left(\sum _j|b_{ij}><b_{ij}|\right)\}-Tr\left[\lambda (R+R^{+})\right]$$

$$Q_i=\sum _j|b_{ij}><a_{ij}|$$

$$S_Q=\sum _i||G_i-Q_i|{|}^2\text{ (for problem 2a)}$$

$$S_Q=\sum _{ij}||G_i{a}_{ij}-b_{ij}|{|}^2\text{ (for problem 2b)}$$

Full text

Next we considered small unitary rotation of entire group which makes the group consistent with a set of additional conditions. Again, iterative procedure based on this correction shows fast convergence. All algorithms developed in this paper were implemented in a Python library which is available as open source software.

1 Introduction

This paper deals with a situation when there is a finite group of unitary matrices (or operators) $G=\{G_{i}\}$ which we need to know, but instead of them only approximate matrices $G_{i}^{\prime}\approx G_{i}$ are available. Besides the obvious disadvantage of not being exact, the set of $G^{\prime}_{i}$ matrices has another serious drawback - it becomes inconsistent with G group multiplication table, is not closed any more with respect to multiplication operation and therefore does not comprise a group. In this paper the following problems are considered and solved:

• •

Problem 1. Perform the multiplication table based group reconstruction, i.e. find the matrices $G_{i}$ that are close to approximate matrices $G^{\prime}_{i}$ but do obey the group multiplication table;

• •

Problem 2a. In addition to problem 1, the matrices $G_{i}$ should be as close as possible to a set of given matrices $Q_{i}$ which are not required to be unitary and non-degenerate;

• •

Problem 2b. In addition to problem 1, the matrices $G_{i}$ should turn a set of vectors $a_{ij}$ into vectors $b_{ij}$, i.e. approximate equality $G_{i}a_{ij}\approx b_{ij}$ should be fulfilled with best possible accuracy.

The solutions of all three problems have been implemented as an open source Python library freely available on GitHub grpcorr .

The example application considered in this paper concerns molecular symmetry when the geometry of the molecule is for various possible reasons distorted from symmetric configuration. These reasons can include thermal movement, presence of external field, spontaneous symmetry break, or be rather computational than physical: not fully converged geometry optimisation. However, the suggested algorithm can find much wider application.

The obvious generalisation of molecular and crystalline symmetry groups are the space groups in higher dimensions. In material science such groups are widely used to explain the structure of quasicrystals and other quasiperiodic solids. After the amazing discovery of a metallic alloy with crystallographically forbidden symmetry in 1984 shechtman it was briefly recognised janssen ; janot that its structure can be viewed as a slice of multidimensional periodic lattice projected into 3D space. Since that unusual properties of quasiperiodic solids continue to attract great attention stein ; delodi ; liu ; hua , including recent discovery of superconductivity in quasicrystals kamiya . An up-to-date review can be found in Ref. qrev . An algorithm capable of producing a group out of set of matrices can potentially facilitate the recognition of quasicrystalline symmetries.

The paper is organised as follows. First, in section 2 the problem 1 is considered. Second, in section 3 problems 2a and 2b are formulated as least-squares problems and common solution is suggested for them. A real example concerning molecular symmetry is provided in section 4 together with efficiency and convergence consideration of suggested method. After that the conclusion follows.

2 Multiplication table based group reconstruction

To start with, let us introduce some notations. The group multiplication table will be denoted $(ij)$, i.e. multiplication of $i$-th and $j$-th group elements produces the element number $(ij)$. The group size is denoted $N$ and the deviation of approximate group matrices from exact ones is $\delta G_{i}=G^{\prime}_{i}-G_{i}$.

Let us introduce also $F_{ij}=G^{\prime}_{i}\cdot G^{\prime}_{j}-G^{\prime}_{(ij)}$ - the matrices characterising the multiplication table violation. Up to the first order of magnitude over $\delta G_{i}$ they can be expressed as $F_{ij}=G_{i}\delta G_{j}+\delta G_{i}G_{j}-\delta G_{(ij)}$. Then we calculate

[TABLE]

and

[TABLE]

where $D_{1}=\sum_{j}\delta G_{j}G_{j}^{-1}$, $D_{2}=\sum_{j}G_{j}^{-1}\delta G_{j}$ and $[A,B]=AB-BA$ stands for the commutator. Taking sum of (1) and (2) allows us to get explicit expression for the deviations:

[TABLE]

By no means can the second term of this equation be expressed through the known values. In fact, it can be safely set to zero and there is a good reason for that. Indeed, consider the unitary rotation of entire group $G^{\prime}_{i}=U^{-1}G_{i}U$ with a unitary matrix $U$ being the same for all $G_{i}$. Such a rotation does not violate the group multiplication table. Up to the first order of magnitude small unitary rotation can be represented as $U=1+R$ with small anti-hermitian $R$. Thus, with the same accuracy

[TABLE]

Here ${}^{\prime}+^{\prime}$ superscript stands for hermitian conjugate. The deviations of the form (4) are not distinguishable for multiplication table based group reconstruction and should be either set to zero or determined from some other sources.

Second term of (3) has exactly this structure. Indeed, since both $G_{i}$ and $G^{\prime}_{i}$ are unitary, $(G^{\prime}_{i})^{+}G^{\prime}_{i}-G_{i}^{+}G_{i}=\delta G_{i}^{+}G_{i}+G_{i}^{+}\delta G_{i}=0$ leading to the expression $\delta G_{i}^{+}=-G_{i}^{-1}\delta G_{i}G_{i}^{-1}$. It is trivial then to show that $[G_{j}^{-1},\delta G_{j}]$ commutators forming second term of (3) are anti-hermitian.

Substituting $G_{i}$ (which are unknown) with $G^{\prime}_{i}$ (which are known) in the expression (3) leads to the error only in the second order of magnitude over $\delta G_{i}$. We have arrived to the final first-order of magnitude expression for the deviations that solves Problem 1:

[TABLE]

The computational procedure based on this expression is simple and straightforward. It is applied iteratively until $G_{i}$ matrices are converged. Each iteration includes:

Re-unitarisation of $G_{i}$ matrices to enforce their unitarity. Re-unitarisation is done according to the formula $G_{i}\rightarrow(G_{i}G_{i}^{+})^{-1/2}G_{i}$; 2. 2.

Construction of multiplication table violation matrices $F_{ij}$; 3. 3.

Calculation of error value characterising the violation of multiplication table

[TABLE]

If the error is smaller than the convergence threshold value $\varepsilon$ then iterations are terminated and $G_{i}$ matrices are considered successfully converged; 4. 4.

Application of formula (5) and subtraction of $\delta G_{i}$ deviations from $G_{i}$ matrices.

The algorithm turns out to be quadratically convergent and converges very fast. The numerical consideration is provided below in section 4.

3 Group correction with simultaneous least squares fit

The above described procedure allows one to find the $G_{i}$ matrices obeying the group multiplication table. However, these matrices are still not uniquely defined. They are defined only up to arbitrary unitary rotation of entire group, and this remaining degree of freedom can be used to solve Problems 2a and 2b. First, we formulate them as least-squares problems:

• •

$\sum_{i}||G_{i}-Q_{i}||^{2}=min$ (problem 2a);

• •

$\sum_{ij}||G_{i}a_{ij}-b_{ij}||^{2}=min$ (problem 2b);

Here $||A||^{2}=\sum_{\alpha\beta}|A_{\alpha\beta}|^{2}=Tr(A^{+}A)$ is the matrix norm squared. Now we consider unitary rotation of the group $G_{i}\rightarrow e^{-R}G_{i}e^{R}$ expressed through small anti-hermitian matrix $R$. Anti-hermiticity condition $R+R^{+}=0$ should be used as a constraint when minimising the sum of squares for problems 2a and 2b. Our goal is to find $R$ matrix minimising the corresponding error function (sum of squares).

Now we are ready to write down the value which is to be minimised in order to solve Problem 2a:

[TABLE]

Here $\lambda$ is the Lagrange multipliers matrix to account for the constraint. In order to obtain linear equation for $R$, this expression has to be expanded up to second order over $R$ and subsequently differentiated:

[TABLE]

Similarly, one can proceed with the differentiation over $R^{+}$:

[TABLE]

Subtraction of equations (8) and (9) allows one to obtain simpler equation, and, what is more important, to get rid of Lagrange multipliers $\lambda$. After that it is safe to substitute $R^{+}\rightarrow-R$ everywhere and get the final equation for R matrix:

[TABLE]

with $H_{i}$ denoting the following structure:

[TABLE]

The matrix equation (3) falls into the class of Generalised Sylvester equations hautus . There is a straightforward but possibly not very efficient solution of such equation. It can be seen as a system of linear equations for elements of matrix R. To re-formulate it in this manner, we need the matrix vectorisation operator $\text{vec}\,R$ that puts all columns of R in the same vector one after another and the Kronecker product operation

[TABLE]

Now (3) is equivalent to the linear equations system $L\cdot\text{vec}\,R=\text{vec}\,P$ with the right-hand side denoted

[TABLE]

and $L$ being the “supermatrix”:

[TABLE]

Here ${}^{\prime}T^{\prime}$ and ${}^{\prime}*^{\prime}$ superscripts stand for the transpose and complex conjugate of matrices, respectively, while $I_{n}$ is the identity matrix of dimension $n$. Formulae (11,12,13) together with the expression $\text{vec}\,R=L^{-1}\cdot\text{vec}\,P$ provide the formal solution of Problem 2a.

However, this solution can be complicated by the presence of zero eigenvalues of $L$ supermatrix. Such zero eigenvalues correspond to eigenmatrices which, being substituted instead of $R$ into left-hand side of (3), turn it to zero. There is always at least one such eigenmatrix with zero eigenvalue - it is the identity matrix. However, it turns out that all such eigenmatrices are hermitian and can be easily ruled out when constructing the solution of (3).

Another serious drawback of the suggested solution is its scalability. Indeed, if the matrices $G_{i}$ have the dimension $n\times n$, then $L$ supermatrix has already $n^{2}\times n^{2}$. Taking into account that commonly used eigensolvers for dense matrices scale cubically, that gives us total scalability of $n^{6}$. For large $n$ it can be computationally very expensive and even prohibitive.

However, a great simplification of this solution can be done if it can supposed that $G_{i}\approx Q_{i}$. Denoting $\delta Q_{i}=Q_{i}-G_{i}$ one can see that right-hand side of (3) depends on these quantities linearly: $P=\sum_{i}[G_{i},\delta Q_{i}^{+}]+[G_{i}^{+},\delta Q_{i}]$. After that substituting $Q_{i}$ with $G_{i}$ on the left-hand side of this equation creates an error only of second order of magnitude over $\delta Q_{i}$.

Then immediately each $H_{i}=2$ because of the unitarity of $G_{i}$ matrices. Using once again the fact that $G_{i}^{+}=G_{i}^{-1}$ and reordering the summation in some terms at the left-hand side of equation (3) gives us:

[TABLE]

Consider the eigenmatrices of left-hand side of this equation. Such an eigenmatrix should obey

[TABLE]

with eigenvalue $\alpha$ which we represent as $\alpha=4(N-\mu)$. Then $\mu$ is an eigenvalue for another matrix transformation: $\sum_{i}G_{i}RG_{i}^{+}=\mu R$. Multiplying by $G_{j}$ on the left side and by $G_{j}^{+}$ on the right side and reordering the summation leads us to the conclusion that $\mu R=\mu G_{j}RG_{j}^{+}$ for any $j$, in other words, that the matrix $\mu R$ must commute with the entire group $G$.

The matrices of the group $G$ themselves comprise some representation of this group. It can be either irreducible or reducible, in that latter case they can be decomposed into the direct sum or irreducible representations: $G_{i}=G_{i}^{(1)}\oplus G_{i}^{(2)}\oplus\dots\oplus G_{i}^{(k)}$. In other words, there is a unitary basis set transformation U which makes all $G_{i}$ matrices block-diagonal with the blocks being the matrices of corresponding irreducible representation. According to the Schur’s lemma any matrix commuting with irreducible representation matrices must be identity matrix multiplied by a scalar.

Thus, we arrived to the conclusion that $\mu R=U^{+}(\mu_{1}I_{n_{1}}\oplus\mu_{2}I_{n_{2}}\oplus\dots\oplus\mu_{k}I_{n_{k}})U$. If this construction is non-zero, then it is necessarily hermitian. For any anti-hermitian matrix $R$ it means that $\mu=0$, and we can get simple explicit expression for $R$:

[TABLE]

This formula, although being approximate as the others in the present paper, is correct up to first order over $\delta Q_{i}$ and is free from the drawbacks of previous supermatrix-based solution.

It remains only to show that Problem 2b can be reduced to Problem 2a. Error function (sum of squares) for Problem 2b is

[TABLE]

If we set

[TABLE]

then second and third terms under the sum over $i$ become identical to those in (7), while the first and the fourth do not, but they make no contribution into equation (3). This means that Problem 2b can be solved identically to Problem 2a with such a choice of $Q_{i}$ matrices.

The algorithm solving these two problems is, again, iterative. Each iteration consists of the following steps:

Re-unitarisation of $G_{i}$ matrices; 2. 2.

Calculation of multiplication table violation error $S_{M}$ given by (6). If this value exceeds the threshold $\varepsilon_{M}$, one step of multiplication table based reconstruction is performed and the convergence flag is set to “False”; 3. 3.

Calculation of $P$ matrix via (12); 4. 4.

Calculation of $R$ by either supermatrix-based approach ($\text{algo}=1$) or by simplified formula (14) ($\text{algo}=2$). The algo parameter is the choice of the user. 5. 5.

Rotation of the group as $G_{i}\rightarrow e^{-R}G_{i}e^{R}$; 6. 6.

Calculation of error value characterising the difference between $G_{i}$ and $Q_{i}$ matrices

[TABLE]

or

[TABLE]

and its change on current iteration $\Delta S_{Q}=S_{Q}(prev)-S_{Q}(current)$; 7. 7.

Comparison of $\Delta S_{Q}$ with the threshold value $\varepsilon_{Q}$ and comparison of $||R||$ with $\varepsilon_{R}$. If one of these values exceeds convergence threshold, the convergence flag is set to “False”; 8. 8.

Checking the convergence flag. If it is “True”, the fit is successfully converged.

4 Application to molecular symmetry

Application of the reported methods to the symmetry of molecules, nanostuctures and crystals has been in fact the main driving force behind this work, although their applicability turned out to be much wider. In this section we consider how the suggested algorithms work for the groups of molecular symmetry. Through this section we use the Shönflies notation for point groups flurry .

First, the multiplication table based group reconstruction (the solution of problem 1) has been tested. Two molecules, namely, C20 fullerene and sulphur hexafluoride were selected for this test. Their symmetry groups are Ih and Oh, respectively. The molecules were constructed as simple list of Cartesian coordinates of each atom, the length unit for the coordinates was Ångström.

Multiplication table for the symmetry group of each molecule has been constructed as follows. All possible permutations of atoms that do not change interatomic distances and bond angles comprise the permutation group isomorphic to the point symmetry group of the molecule. This permutations group was used to build the multiplication table. After that geometries of both molecules were distorted by adding random displacements to each coordinate. Each displacement value for each coordinate ranged uniformly from -0.5 to 0.5 Å. Next, the group matrices of initial approximation were constructed via formula (15) and re-unitarised. Finally, the multiplication table based group reconstruction has been applied. The driver routine implemented in the repository grpcorr is called multab_group_correction.

The results are shown on fig. 1. It requires only 4 iterations to achieve the convergence of $10^{-12}$ for the multiplication table violation error $S_{M}$. Note also that the convergence becomes faster with each next iteration due to the quadratucally convergent character of the suggested algorithm.

The insets on both figs. 1a and 1b show the randomly distorted molecule before reconstruction and symmetrised molecule after the group reconstruction. The rotation axis of the symmetry group are depicted with arrows. The highest order axis (5-th for C20 and 4-th for SF6) are shown black while the others have different colours. Observe the bundles of axis before the reconstruction turns them into just one axes. Due to inconsistency of approximate group some rotations which should have exactly the same axes have close but different ones. Reconstruction algorithm successfully corrects both rotation axis and angles.

Next, the least-squares fit algorithm has been tested. As an example fig. 2 shows the ethane molecule. In the beginning the molecule had its full symmetry D3d which was correctly determined similarly to previous example. Then C2H6 molecule was rotated around a random axis by random angle (on fig. 2 it was 86.47*∘*). After that the algorithm solving Problem 2b has been applied (the corresponding driver routine in the repository grpcorr is called abfit_group_correction). The convergence on fig. 2 is characterised by previously defined quantities $S_{Q}$ calculated according to (17), its change and the rotation norm $||R||$. As one can see, all three quantities decay exponentially during the fit and after it is finished the symmetry group axes are oriented properly with respect to the molecule.

Finally, we have tested the ability of suggested algorithms to manage both types of error simultaneously, i.e. to correct the group multiplication table and at the same time find the unitary rotation. This is illustrated by fig. 3. CH4 molecule has been randomly deformed (each coordinate received random displacement from -0.5 to 0.5 Å). Then initial group matrices were constructed by formula (15) and re-unitarised. After that the molecule was randomly rotated (by 56.93*∘* on fig. 3). When all this preparations were finished, the algorithm for Problem 2b has been applied.

The convergence was fastest for $S_{M}$ value similar to previous examples. After only 4 iterations it became virtually zero. The values of $\Delta S_{Q}$ and $||R||$ showed exponential decay and converged to threshold value during 24 iterations. In this example the value of $S_{Q}$ does not converge to zero because the molecule is not only rotated but also deformed. Instead it converges to its minimal value (8.316 for the example shown on fig. 3).

It is interesting to note that algorithm for Problem 2b always finds the correct rotation only when the rotation angle is less then 90*∘*. For larger angles the convergence is never achieved. As the test calculations described in this section were repeated several hundreds of times this can be stated with confidence. Another interesting observation is that there is almost no difference in the convergence between the “supermatrix” approach and simplified formula (14) for the problems 2a and 2b.

5 Conclusion

In this paper we have considered a finite unitary matrix group with exact elements being unknown and only approximate elements available. The first problem of such a group is the inconsistency with its own multiplication table which we called “Problem 1”. Another distinct inconsistency type arises when it is known how the group elements should act on certain set of vectors, but they do not act as they should. We formulated this as “Problem 2a” and “Problem 2b”.

The solutions of all three considered problems were found analytically as a first order correction to the group elements. Those solutions can be used in iterative procedure until the group matrices are converged, and such iterative procedures were implemented in the GitHub repository grpcorr . The present Python implementation does not pretend to be highly computationally efficient. It was programmed rather to test the suggested methods and explore their applicability. However, as it is published as open-source software, any researcher can feel free to use the code or to rewrite it in more efficient manner.

The numerical test performed for the point symmetry groups had shown very good convergence of suggested methods. The multiplication table based reconstruction has quadratic convergence, while for the least-squares fit of unitary rotation it is linear. Although initially we considered the situation when the deviations of approximate group matrices from exact ones are small, the suggested algorithms seem to successfully handle the situations when they are not.