## Worked example: Symmetry reduction


```python
import numpy as np
import sympy as sp
import scipy.linalg as la
```

### Some routines for creating rotation matrices and transpositions


```python
def get_rotation_matrix(phi, axis):
    cp = np.cos(phi)
    sp = np.sin(phi)
    r1, r2, r3 = axis
    Rm = np.array(
        [
            [
                r1 ** 2 * (1 - cp) + cp,
                r1 * r2 * (1 - cp) - r3 * sp,
                r1 * r3 * (1 - cp) + r2 * sp,
            ],
            [
                r1 * r2 * (1 - cp) + r3 * sp,
                r2 ** 2 * (1 - cp) + cp,
                r2 * r3 * (1 - cp) - r1 * sp,
            ],
            [
                r3 * r1 * (1 - cp) - r2 * sp,
                r2 * r3 * (1 - cp) + r1 * sp,
                r3 ** 2 * (1 - cp) + cp,
            ],
        ]
    )
    # Clean spurious terms
    for ii, jj in np.ndindex(Rm.shape):
        if abs(Rm[ii, jj]) < 1e-10:
            Rm[ii, jj] = 0
    sp = np.sin(phi)
    return Rm


def get_transposition(nd, verbose=0):
    # Make a basis:
    basis_nd = []
    for ii in range(nd):
        for jj in range(nd):
            emp = np.zeros([nd, nd], dtype="int8")
            bb = emp
            bb[ii, jj] = 1
            basis_nd.append(bb)
    #
    # Build transpose and return:
    transp = np.array([bb.T.flatten() for bb in basis_nd])

    return transp


# pretty print matrix
def pprint(M, maxl=80):
    M = np.atleast_2d(M)
    dim1, dim2 = M.shape
    lengths = [len(str(el)) for el in M.flatten()]
    maxl = min(maxl, max(lengths))
    for ii in range(dim1):
        print(
            (f" | " + " ,  ".join([f"{str(el)[:maxl]:{maxl}s}" for el in M[ii]]) + " |")
        )
```

### Start: Define a matrix $A$
* Define a general 3x3 matrix $A$ and its 9x1 vectorized representation
$$ \boldsymbol{a} = {\textbf{vec}} (A)$$

* $A$ and $\boldsymbol a$ can, e.g., represent a physical tensor (conductivity, stress, ...)


```python
ndim = 3
sym = sp.Symbol
A = sp.ones(ndim, ndim)
for ii in range(ndim):
    for jj in range(ndim):
        A[ndim * ii + jj] = sym(f"a_{ii}{jj}")
A = np.array(A)
a = A.flatten()
print("Matrix A:")
pprint(A)
print("vec(A):")
pprint(a)
```

    Matrix A:
     | a_00 ,  a_01 ,  a_02 |
     | a_10 ,  a_11 ,  a_12 |
     | a_20 ,  a_21 ,  a_22 |
    vec(A):
     | a_00 ,  a_01 ,  a_02 ,  a_10 ,  a_11 ,  a_12 ,  a_20 ,  a_21 ,  a_22 |


### Case 1: Cubic system
* Define 3 rotation matrices $M_i$ that implement a 90° rotation about $x$, $y$, and $z$ axis
* Construct the rotation matrices in the flattened representatin by Roth's Relationship, i.e.
$$ M_\text{flat} = M \otimes M $$
* Add transposition (not necessary in cubic case)


```python
r1 = get_rotation_matrix(np.pi / 2, [1, 0, 0])
r2 = get_rotation_matrix(np.pi / 2, [0, 1, 0])
r3 = get_rotation_matrix(np.pi / 2, [0, 0, 1])
pprint(r1)
```

     | 1.0  ,  0.0  ,  0.0  |
     | 0.0  ,  0.0  ,  -1.0 |
     | 0.0  ,  1.0  ,  0.0  |


#### Construct big matrices implementing rotations (+transposition) of the vectorized tensor


```python
R1 = np.kron(r1, r1)
R2 = np.kron(r2, r2)
R3 = np.kron(r3, r3)
T  = get_transposition(ndim)
```

#### Now sum up the matrices we want to invariant under, i.e.
$$ \sum_i (\mathbf 1 - M_i)~a = 0$$


```python
id = np.eye(len(a))
inv = 4*id - R1 - R2 - R3 - T
```

#### Consctruct nullspace by SVD
__[scipy-cookbook.readthedocs.io/items/RankNullspace.html](http://scipy-cookbook.readthedocs.io/items/RankNullspace.html)__


```python
u, s, vh = la.svd(inv, lapack_driver="gesdd")
rank = (s > 1e-12).sum()

print(f"Initial Dimension:                              {len(a)}")
print(f"Rank of Nullspace (= No. irred. elements):      {len(a) - rank}")
```

    Initial Dimension:                              9
    Rank of Nullspace (= No. irred. elements):      1


#### Construct matrices translating between full and reduced representation

* $S$ reconstructs the full representation $a$ from a given reduced $\tilde a$. One can think of $\tilde a$ as being the components of $a$ in the basis given by the vectors in $S$:
$$ a = S \, \tilde{a}$$

* $S^+$ (pseudo-inverse, often just transpose) projects a given $a$ onto the irreducible components $\tilde a$ 
$$\tilde{a} = S^+ a$$


```python
S = vh[rank:].T
Sp = S.T
```

#### Build projectors
* $P = S^\vphantom{+} S^+$

* $\boldsymbol{1}_\text{irred} = S^+ S^\phantom{+}$


```python
P = S @ Sp
id_irred = Sp @ S
print("Projector onto invariant space")
pprint(P, 4)
print("Identity within invariant space")
pprint(id_irred)
```

    Projector onto invariant space
     | 0.33 ,  0.0  ,  1.77 ,  -6.9 ,  0.33 ,  1.01 ,  1.29 ,  -5.2 ,  0.33 |
     | 0.0  ,  0.0  ,  0.0  ,  0.0  ,  0.0  ,  0.0  ,  0.0  ,  0.0  ,  0.0  |
     | 1.77 ,  0.0  ,  9.49 ,  -3.6 ,  1.77 ,  5.44 ,  6.91 ,  -2.8 ,  1.77 |
     | -6.9 ,  0.0  ,  -3.6 ,  1.43 ,  -6.9 ,  -2.1 ,  -2.6 ,  1.09 ,  -6.9 |
     | 0.33 ,  0.0  ,  1.77 ,  -6.9 ,  0.33 ,  1.01 ,  1.29 ,  -5.2 ,  0.33 |
     | 1.01 ,  0.0  ,  5.44 ,  -2.1 ,  1.01 ,  3.11 ,  3.96 ,  -1.6 ,  1.01 |
     | 1.29 ,  0.0  ,  6.91 ,  -2.6 ,  1.29 ,  3.96 ,  5.04 ,  -2.0 ,  1.29 |
     | -5.2 ,  0.0  ,  -2.8 ,  1.09 ,  -5.2 ,  -1.6 ,  -2.0 ,  8.28 ,  -5.2 |
     | 0.33 ,  0.0  ,  1.77 ,  -6.9 ,  0.33 ,  1.01 ,  1.29 ,  -5.2 ,  0.33 |
    Identity within invariant space
     | 1.0000000000000002 |


#### Symmetrize the tensor with the projector obtained
$$ \text{sym} (a) = P a = S^\vphantom{+} S^+ a = S \, \tilde a $$


```python
aT = np.dot(Sp, a)
aS = np.dot(S, aT) # = np.dot(P, a)
```

#### How does the matrix $A$ now look like?
$$A = \text{unvec} \left( \text{sym} (a) \right)$$


```python
As = aS.reshape(3,3)

print('1/3*')
pprint(3*As)
```

    1/3*
     | 1.0*a_00 + 5.33729362479519e-33*a_02 - 2.07749100017982e-35*a_10 + 1.0*a_11 + 3. ,  0                                                                                ,  5.33729362479519e-33*a_00 + 2.84867032372794e-65*a_02 - 1.10881794708291e-67*a_1 |
     | -2.07749100017982e-35*a_00 - 1.10881794708291e-67*a_02 + 4.31596885582815e-70*a_ ,  1.0*a_00 + 5.33729362479519e-33*a_02 - 2.07749100017982e-35*a_10 + 1.0*a_11 + 3. ,  3.05816280424746e-34*a_00 + 1.63223128386957e-66*a_02 - 6.35330570290877e-69*a_1 |
     | 3.8891592043194e-34*a_00 + 2.07575846270274e-66*a_02 - 8.07969324524005e-69*a_10 ,  -1.57643859172956e-33*a_00 - 8.41391564551927e-66*a_02 + 3.2750369866543e-68*a_1 ,  1.0*a_00 + 5.33729362479519e-33*a_02 - 2.07749100017982e-35*a_10 + 1.0*a_11 + 3. |


$= \frac{1}{3} \text{Tr A}$ 

## How about hexagonal?
Start with defining three-fold rotations in $x,y$ plane, i.e., about $z$ axis


```python
r1 = get_rotation_matrix(np.pi / 6, [0, 0, 1])
r2 = get_rotation_matrix(np.pi / 3, [0, 0, 1])
```

Construct the matrices in the flatten represenation as before


```python
R1 = np.kron(r1, r1)
R2 = np.kron(r2, r2)
T = get_transposition(ndim)
```


```python
id = np.eye(len(a))
inv = 3*id - R1 - R2 - T
```


```python
# Consctruct nullspace
u, s, vh = la.svd(inv, lapack_driver="gesdd")
rank = (s > 1e-12).sum()

print(f"Dimension:          {len(a)}")
print(f"Rank of Nullspace:  {len(a) - rank}")
```

    Dimension:          9
    Rank of Nullspace:  2



```python
# Nullspace:
S = vh[rank:].T
# clean S
for ii, jj in np.ndindex(S.shape):
    if abs(S[ii, jj]) < 1e-10:
        S[ii, jj] = 0
Sp = S.T
```

How do the projectors look like?


```python
print(S@Sp)
#print(Sp)
#print(S@Sp)
print(Sp@S)
```

    [[0.5 0.  0.  0.  0.5 0.  0.  0.  0. ]
     [0.  0.  0.  0.  0.  0.  0.  0.  0. ]
     [0.  0.  0.  0.  0.  0.  0.  0.  0. ]
     [0.  0.  0.  0.  0.  0.  0.  0.  0. ]
     [0.5 0.  0.  0.  0.5 0.  0.  0.  0. ]
     [0.  0.  0.  0.  0.  0.  0.  0.  0. ]
     [0.  0.  0.  0.  0.  0.  0.  0.  0. ]
     [0.  0.  0.  0.  0.  0.  0.  0.  0. ]
     [0.  0.  0.  0.  0.  0.  0.  0.  1. ]]
    [[1. 0.]
     [0. 1.]]



```python
# Projector
P = S@Sp
```


```python
# Symmetrize a
aS = np.dot(P, a)
```


```python
# Restore shape
As = aS.reshape(3,3)
```


```python
pprint(As)
```

     | 0.5*a_00 + 0.5*a_11 ,  0                   ,  0                   |
     | 0                   ,  0.5*a_00 + 0.5*a_11 ,  0                   |
     | 0                   ,  0                   ,  1.0*a_22            |



```python

```

7c7fb3f348760336303059077dc592f38b4ae7d1
ipynb
Jupyter Notebook
mathematics/nullspace/coffee_talk18.ipynb
flokno/python_recipes
a09da65528ce3a2f1fd884aea361275b9aab0c15
0BSD
2020-05-06T14:38:16.000Z
2022-01-05T09:43:04.000Z
mathematics/nullspace/coffee_talk18.ipynb
flokno/python_recipes
a09da65528ce3a2f1fd884aea361275b9aab0c15
0BSD
mathematics/nullspace/coffee_talk18.ipynb
flokno/python_recipes
a09da65528ce3a2f1fd884aea361275b9aab0c15
0BSD
Qwen/Qwen-72B
1. YES
2. YES