Offline part

Dec 13, 2024 · 17 min read

OFFLINE PART (SVD/POD/Greedy/rSVD…)

Script by Elise Grosjean (elise.grosjean@ensta-paris.fr)

import sys
!{sys.executable} -m pip install numpy
!{sys.executable} -m pip install matplotlib
!{sys.executable} -m pip install scikit-fem

[notice] A new release of pip is available: 24.1.1 -> 25.0.1
[notice] To update, run: python3.11 -m pip install --upgrade pip

Let us consider a parameterized problem we want to solve for different parameter values in a set $\mathcal{G}$.

In general, Reduced Basis Methods (RBM) are decomposed in two stages:

  • one offline that generates a smaller subspace $X_h^N$ of dimension $N$,
  • and one online that find an approximation for a new parameter $\mu \in \mathcal{G}$ onto that subspace.

The offline stage is usually very costly but only executed once. In what follows, we review different algorithms that we can use to compute a reduced basis of $X_h^N$ and compare their computational costs.

# import packages
import skfem  # for Finite Element Method
import numpy as np
import matplotlib.pyplot as plt
from sklearn.utils.extmath import randomized_svd
import time

1) The 2D-driven cavity problem:

We are going to use in this example the famous 2D-driven cavity problem with the Finite Element Method (FEM), which consists in solving on a unit square (denoted $\Omega$) the following equations: \begin{align} &-\nu \Delta u + (u \cdot \nabla) u + \nabla p =0, \textrm{ in } \Omega,\ & \nabla. u=0, \textrm{ in } \Omega,\ & (u_1,u_2)=(1,0), \textrm{ on } \Omega_{up}:=\partial \Omega \cap {y=1},\ & (u_1,u_2)=(0,0), \textrm{ on } \partial \Omega \backslash \Omega_{up}, \end{align} where $u=(u_1,u_2) \in V:=H^1_{d,0}(\Omega)^2=\{u \in H^1(\Omega)^2, \gamma_{\partial \Omega \backslash \Omega_{up}} u=0, \ \gamma_{up} u =(1,0) \}$ ($\gamma$ stands for the trace operator) represents the velocity of the incompressible fluid, $ p \in L^2(\Omega)$ its pressure, and $\nu=\frac{1}{Re}$ where $Re$ is the Reynolds parameter. Here, the Reynolds number is our parameter of interest ($\mu=Re$). For the nonlinearity we adopt a fixed-point iteration scheme, and after multiplying by test functions $q$ and $v$ (resp. for pressure and velocity), which in variational form reads: \begin{equation} \nu (\nabla uk, \nabla v) + ((u{k-1} \cdot \nabla) uk,v) - (pk, \nabla \cdot v) - (q, \nabla \cdot uk) + 10{-10} (p^k, q) =0, \textrm{in } \Omega, \end{equation} where $u^{k-1}$ is the previous step solution, and we iterate until a threshold is reached (until $\|u^{k}-u^{k-1}\| < \varepsilon $). Here, with the term $ 10^{-10} (p^k, q) $, we impose the average of the pressure $\int_{\Omega} p^k $ to be equal to $0$. For more details on how we derive this formulation, visit the link : https://github.com/grosjean1/navierStokes (FEM.pdf).

We employ Taylor-Hood elements to get a proper solution (e.g. P2-P1 for the tuple velocity-pressure) and obtain the system $\mathbf{K} \mathbf{x} =\mathbf{f}$ to solve where $\mathbf{K}= \begin{pmatrix} \mathbf{A} & -\mathbf{B}^T\\ -\mathbf{B} & 10^{-10} \mathbf{C} \end{pmatrix}$, $\mathbf{x}$ stands for the tuple velocity-pressure $(u_1^k,u_2^k,p^k)$, and where the assembled matrix $\mathbf{A}$ corresponds to the bilinear part $ \nu (\nabla u^k, \nabla v) + ((u^{k-1} \cdot \nabla) u^k),v) $, the matrix $ B$ to the bilinear part $( p^k ,\nabla \cdot v)$ and $\mathbf{C}$ is the mass matrix applied to the pressure variable ($(\cdot,\cdot)$ represents either the $L^2$ inner-product onto the velocity space or onto the pressure space).

The dirichlet boundary conditions are imposed with a penalization method, called with scikit-FEM by the line:

solve(*condense(K,f, x=uvp, D=D)),

where x=uvp gives the values at the boundaries of the velocity and D refers to the boundary decomposition.

# First we define a mesh for the unit square
mesh= skfem.MeshTri().refined(5).with_boundaries(                                                                
        {                                                                                                                                
            "left": lambda x: x[0] == 0,                                                                                            
            "right": lambda x: x[0] == 1,            
            "down": lambda x: x[1] == 0,                                                                                            
            "up": lambda x: x[1] == 1,     
            
        }                                                                                                                               
)

print(mesh)
mesh

<skfem MeshTri1 object>
  Number of elements: 2048
  Number of vertices: 1089
  Number of nodes: 1089
  Named boundaries [# facets]: left [32], bottom [32], right [32], top [32], down [32], up [32]

svg

"""First step."""

# Assembling matrices

from skfem.assembly import BilinearForm, LinearForm
from skfem.helpers import grad, dot,div,ddot

@BilinearForm
def laplace(u, v, _):
    print("gradula",np.shape(grad(u)))
    return dot(grad(u), grad(v))

@BilinearForm
def vector_laplace(u, v, w):
    # same as laplace but for u,v vectors 
    return w.nu*ddot(grad(u), grad(v))

@BilinearForm
def mass(u, v, _):
    return u * v

@BilinearForm
def divu(u, v, w):
    return div(u) * v
divergence = divu
 
    #z=w.up*u.grad[0] + w.up * u.grad[1] #(u^(k-1). grad u^k). v
@BilinearForm
def nonlinearterm(u, v, w):
    up_x, up_y = w.up 
    v_x,v_y=v[0],v[1]
    
    gradu=grad(u)
    gradxu_x=gradu[0,0] 
    gradxu_y=gradu[0,1]
    gradyu_x=gradu[1,0]
    gradyu_y=gradu[1,1]
    
    return (up_x*gradxu_x+up_y*gradxu_y)*v_x + (up_x*gradyu_x+up_y*gradyu_y)*v_y

In the 2D-driven cavity problem, the parameter of interest is the Reynolds number. In the code below, the function SolveCavityProblem(Re,Mesh) takes as parameters the Reynolds number and a mesh, and returns the associated solution (velocity-pressure).

""" Finally we solve the problem """
from skfem import *
# Compute solution of the 2d-driven cavity problem
element = {'u': ElementVector(ElementTriP2()),
           'p': ElementTriP1()} #[u,v,p]= velocity-pressure with Taylor-Hood P2-P1 FEM elements


def SolveCavityProblem(Re,Mesh):
    # Re is the reynolds parameter #
    Nu=1./Re 
    print('Parameter nu:',Nu)

    basis = {variable: Basis(Mesh, e, intorder=4)
         for variable, e in element.items()} # FEM space
    up=basis['u'].zeros() #initialization for non-linear term previous solution
    
    B =divu.assemble( basis['u'], basis['p']) # B.T=p*div(v) and B=q *div(u) with q=pressure test function
    C = mass.assemble(basis['p'],basis['p'] ) # 1^-10 * p*q impose pressure average equal to 0.
    A1 =vector_laplace.assemble(basis['u'],nu=Nu) # nu* grad(u) * grad(v) with u=[u1,u2] and v is test function v=[v1,v2]
    #A2=nonlinearterm.assemble(basis['u'],up=basis['u'].interpolate(up)) #(u^(k-1). grad u^k). v


    # global matrix Stokes
    K = bmat([[A1, -B.T],
          [-B, 1e-10 * C]], 'csr')

    def profil_up(x):
        """return u=(0,1) at the up boundary and (0,0) otherwise """
        return np.stack([x[1]==1, np.zeros_like(x[0])])
    
    up_basis = FacetBasis(Mesh, element['u'], facets=Mesh.boundaries['up'])
    all_basis = FacetBasis(Mesh, element['u'])

    uvp_boundary = np.hstack((
        all_basis.project(profil_up), # can be replaced by all_basis.project(profil_up), will be 0 by default for the other boundaries
        basis['p'].zeros(),
    ))

    D = basis['u'].get_dofs(['up', 'down', 'right','left']) # since we put D on basis['u'], there will be no penalization on pressure
    f = np.concatenate([basis['u'].zeros(),basis['p'].zeros()]) # =0

    #######################
    ### SOLVING PROBLEM ###
    #######################
    uvp = solve(*condense(K,f, x=uvp_boundary, D=D))
    velocity, pressure = np.split(uvp, K.blocks)

    for i in range(15): #fixed-point iterations
        up=velocity.copy() # updating velocity
        A2=nonlinearterm.assemble(basis['u'],up=basis['u'].interpolate(up))

     
        # global matrix assembling update
        K = bmat([[A1+A2, -B.T],
          [-B, 1e-10 * C]], 'csr')

        uvp = solve(*condense(K,f, x=uvp_boundary, D=D))
     
        velocity, pressure = np.split(uvp, K.blocks)

    return velocity,pressure
    

""" We may visualize the solutions """

from skfem.visuals.matplotlib import plot, draw, savefig

velocity,pressure=SolveCavityProblem(150,mesh) 

def visualize_pressure(Mesh,Pressure):
    ## Function that allows us to visualize the pressure
    #from skfem.visuals.matplotlib import plot
    return plot(Mesh, Pressure, shading='gouraud', colorbar=True)
    
visualize_pressure(mesh,pressure).show()

def visualize_velocity(Mesh,Velocity):
    ## Function that allows us to visualize the velocity
    fig, ax = plt.subplots()

    # Calculate the magnitude of velocity
    basis = {variable: Basis(Mesh, e, intorder=4)
            for variable, e in element.items()} # FEM space
    velocity1 = Velocity[basis['u'].nodal_dofs].copy()
    magnitude = np.sqrt(velocity1[0]**2 + velocity1[1]**2)

    # Display the velocity fields with colors for magnitude
    quiver = ax.quiver(*mesh.p, *velocity1, magnitude, angles='xy', cmap='jet')

    # colorbar
    cbar = plt.colorbar(quiver, ax=ax)

    # Limits of axis for the mesh
    ax.set_xlim([Mesh.p[0].min(), Mesh.p[0].max()])
    ax.set_ylim([Mesh.p[1].min(), Mesh.p[1].max()])
    cbar.set_label('Velocity magnitude')
    # show results
    plt.show()

visualize_velocity(mesh,velocity)

Parameter nu: 0.006666666666666667

png

png

2) The Reduced Basis (RB) construction:

For all the methods listed below, we need a fine mesh and the generation of snapshots: We decompose the snapshots by one average over the parameters and by one fluctuation part. Then, the modes are estimated with the fluctuations.

## FINE MESH

FineMesh = skfem.MeshTri().refined(4).with_boundaries(                                                                
        {                                                                                                                                
            "left": lambda x: x[0] == 0,                                                                                            
            "right": lambda x: x[0] == 1,            
            "down": lambda x: x[1] == 0,                                                                                            
            "up": lambda x: x[1] == 1,     
            
        }                                                                                                                               
)
FineBasis = {variable: Basis(FineMesh, e, intorder=4)
            for variable, e in element.items()} # FEM space

NumberOfNodesFineMesh = FineMesh.p.shape[1]
print("number of fine nodes: ",NumberOfNodesFineMesh)
#num_dofs_uFineMesh = basis['u'].doflocs.shape[1] # or np.shape(velocity)[0] for DOFs

number of fine nodes:  289

""" Snapshots generation"""

NumberOfSnapshots=10 #Training set

#####  Create fluctuation parts for the snapshots #####
Re=1. #first Reynolds number
AveragedSnapshots=FineBasis["u"].zeros()

FineSnapshots=[]
for i in range(NumberOfSnapshots):
    velocity,pressure=SolveCavityProblem(Re,FineMesh)
    AveragedSnapshots+=velocity
    FineSnapshots.append(velocity)
    Re+=45

AveragedSnapshots/=NumberOfSnapshots

for i in range(NumberOfSnapshots):
    FineSnapshots[i]-=AveragedSnapshots
    #visualize_velocity(FineMesh,FineSnapshots[i])

@BilinearForm
def massVelocity(u, v, _):
    return u[0] * v[0]+u[1] * v[1]
    
L2=massVelocity.assemble(FineBasis["u"])

Parameter nu: 1.0
Parameter nu: 0.021739130434782608
Parameter nu: 0.01098901098901099
Parameter nu: 0.007352941176470588
Parameter nu: 0.0055248618784530384
Parameter nu: 0.004424778761061947
Parameter nu: 0.0036900369003690036
Parameter nu: 0.0031645569620253164
Parameter nu: 0.002770083102493075
Parameter nu: 0.0024630541871921183

1. Snapshot Proper Orthogonal Decomposition (Snapshot-POD)

One of the most popular algorithm is called “Snapshot-POD”. See details in the paragraph “OFFLINE STAGE” on https://reducedbasis.github.io/docs/pod/.

### SNAPSHOT-POD 
start_time = time.time()

NumberOfModes=6 #tol
print("number of modes: ",NumberOfModes)


# We first compute the correlation matrix C_ij = (u_i,u_j)
CorrelationMatrix = np.zeros((NumberOfSnapshots, NumberOfSnapshots))
for i, snapshot1 in enumerate(FineSnapshots):
    MatVecProduct = L2.dot(snapshot1)
    for j, snapshot2 in enumerate(FineSnapshots):
        if i >= j:
            CorrelationMatrix[i, j] = np.dot(MatVecProduct, snapshot2)
            CorrelationMatrix[j, i] = CorrelationMatrix[i, j]

#print("CorrelationMatrix",CorrelationMatrix)
# Then, we compute the eigenvalues/eigenvectors of C 
EigenValues, EigenVectors = np.linalg.eigh(CorrelationMatrix, UPLO="L") # C eigenVectors=eigenValues eigenVectors
idx = EigenValues.argsort()[::-1] # sort the eigenvalues

TotEigenValues = EigenValues[idx]
TotEigenVectors = EigenVectors[:, idx]

EigenValues=TotEigenValues[0:NumberOfModes] # retrieve only the first eigenvalues and associated eigenvectors
EigenVectors=TotEigenVectors[:,0:NumberOfModes]

print("eigenvalues: ",EigenValues)

RIC=1-np.sum(EigenValues)/np.sum(TotEigenValues) #must be close to 0
print("Relativ Information Content (must be close to 0):",RIC)
#RIClist=[]
#for i in range(len(EigenValues)):
#    RIClist.append(1-np.sum(EigenValues[:i])/np.sum(TotEigenValues)) #must be close to 0

#project the snapshots onto the dominant eigenvectors and orthonormalize to retrieve the reduced basis
ChangeOfBasisMatrix = np.zeros((NumberOfModes,NumberOfSnapshots))
for j in range(NumberOfModes):
    ChangeOfBasisMatrix[j,:] = EigenVectors[:,j]/np.sqrt(EigenValues[j])

ReducedBasis = np.dot(ChangeOfBasisMatrix,FineSnapshots)
end_time = time.time()
execution_time = end_time - start_time
print(f"Temps d'exécution : {execution_time:.6f} secondes")
# orthogonality test:
#for i in range(NumberOfModes):
#    MatVecProduct = L2.dot(ReducedBasis[i])
#    for j in range(NumberOfModes):    
#        test = np.dot(MatVecProduct, ReducedBasis[j])
#        print("i,j (",i,j,"): ",test)

number of modes:  6
eigenvalues:  [2.26782714e-02 2.92187702e-03 2.23860208e-04 1.86248477e-05
 2.17765951e-06 1.26997357e-07]
Relativ Information Content (must be close to 0): 1.9637265658012382e-07
Temps d'exécution : 0.020321 secondes

print("-----------------------------------")
print("  Reduced basis accuracy           ")
print("-----------------------------------")
### Offline Errors (H^1_0 and L2 norm)

L2compressionErrors=[]
H10compressionErrors=[]
H10=vector_laplace.assemble(FineBasis['u'],nu=1)
for snap in FineSnapshots:
    ExactSolution =snap
    CompressedSolutionU= ExactSolution@(L2@ReducedBasis.transpose()) #snapshot compression
    ReconstructedCompressedSolution = np.dot(CompressedSolutionU, ReducedBasis)
    
    norml2ExactSolution=np.sqrt(ExactSolution@(L2@ExactSolution))
    normh1ExactSolution=np.sqrt(ExactSolution@(H10@ExactSolution))
    t=np.abs(ReconstructedCompressedSolution-ExactSolution)
    
    if norml2ExactSolution !=0 and normh1ExactSolution != 0:
        relL2Error=np.sqrt(t@L2@t)/norml2ExactSolution
        relH10Error=np.sqrt(t@H10@t)/normh1ExactSolution
    else:
        relL2Error = np.sqrt(t@L2@t) #np.linalg.norm(ReconstructedCompressedSolution-ExactSolution)
        relH10Error=np.sqrt(t@H10@t)
    L2compressionErrors.append(relL2Error)
    H10compressionErrors.append(relH10Error)

print("L2 compression error =", L2compressionErrors)
print("H1 compression error =", H10compressionErrors)

-----------------------------------
  Reduced basis accuracy           
-----------------------------------
L2 compression error = [1.3815241858207787e-05, 0.00010641993810965995, 0.0004383338831715788, 0.0012111542561164227, 0.001562010925688369, 0.000402636712347126, 0.0008616253925892823, 0.0001517154231303115, 0.000608298618732441, 0.0002466646884607357]
H1 compression error = [1.939956772653135e-05, 0.00014998512464928755, 0.0006436539912057455, 0.0018634474193047483, 0.0024692938629231388, 0.0006456797464323044, 0.0012552738118769562, 0.00023166001504574293, 0.0008375179988110922, 0.0003313113363450415]

cost

If the matrix L2 is dense, then the cost of the correlation matrix computation is in $O(\mathcal{N}_h^2\ Ntrain^2)$ flops ($O(\mathcal{N}_h\ Ntrain^2)$ flops for a sparse L2 matrix), where $\mathcal{N}_h$ is the number of degrees of freedom of one snapshot and $Ntrain$ is the number of snapshots (L2.dot(snapshot1) is done in $O(\mathcal{N}_h^2$ within two imbricated forloops on the number of snapshots). Then, the eigenproblem of the correlation matrix $C$ is solved in $O (\mathcal{Ntrain^3})$ flops.

2. Greedy algorithm

A greedy algorithm is a procedure which aims at approximating each element of a compact set (here the manifold $S_h=Span\{u_h^0, \dots,u_h^n, \}$) in a Hilbert space $V_h$ by a subspace of properly selected elements of this compact set. The greedy procedure is a fast way to compute the modes by choosing some suitable parameters with respect to a criterion. At each Greedy iteration, the new RB parameter is chosen such that the corresponding snapshot is the worse approached by the previous reduced space. The projection of a solution onto the reduced basis is given by $P^n(u)=\underset{i=1}{\overset{n}{\sum}} (u,\Phi_i)\Phi_i$.

  • During the greedy algorithm, we first set $\Phi_0:=\frac{u_h^0}{\|u_h^0\|}$ (or any snapshots in the training set).

  • Then, we look for another snapshot $u_h^k$ in the training set that maximizes the error between itself and its projection onto the reduced basis, built from the first snapshot: $\frac{\|u_h^k - P^0(u_h^k) \|}{\|u_h^k\|}$, and we set $\Phi_1= \frac{u_h^k - P^0(u_h^k)}{\|u_h^k - P^0(u_h^k)\|}$.

  • We iterate until reaching a tolerance threshold on this error maximization, ensuring that any subsequent snapshot projection will result in a smaller error.

start_time = time.time()
ReducedBasis=np.zeros((NumberOfModes,np.shape(velocity)[0]))
SnapshotNorm=np.sqrt(np.dot(FineSnapshots[0],L2.dot(FineSnapshots[0])))

ReducedBasis[0,:]=FineSnapshots[0]/SnapshotNorm #first (random) mode
ListReducedOrderBasis=[ReducedBasis[0,:]]
ListeIndex=[0] #first snapshot 
    
MatrixBasisProduct=[L2.dot(ListReducedOrderBasis[0])]
    
for n in range(1,NumberOfModes):
    
    TestVector=dict() # dictionnary: to find the worst error between the snapshots and their projection on previous RB 
    for j in range(NumberOfSnapshots): 
        if not (j in ListeIndex): #if index not already in the RB 
        
            #w= u_j - sum_k (u_j,Phi_k) Phi_k
            w=FineSnapshots[j]-sum((bk*np.dot(MatrixBasisProduct[k],FineSnapshots[j]) for k,bk in enumerate(ListReducedOrderBasis)))#,axis=0)#potential vector to add in the reduced basis
                
            
            NormL2W=np.sqrt(np.dot((L2.dot(w)),w))
            if np.sqrt(np.dot(FineSnapshots[j],L2.dot(FineSnapshots[j]))) > 1e-10:
                GreedyMaximumTest=NormL2W/np.sqrt(np.dot(FineSnapshots[j],L2.dot(FineSnapshots[j]))) #we seek the max
            else:
                GreedyMaximumTest=NormL2W #we seek the max
            
            TestVector[j]=[GreedyMaximumTest,w,NormL2W]
      
    Index=max(TestVector, key = lambda k: TestVector[k][0]) #index of the snapshot used
    ListeIndex.append(Index) #adding the corresponding j in the list
    
    ListReducedOrderBasis.append(TestVector[Index][1]/TestVector[Index][2])  #orthonormalization in L2
    MatrixBasisProduct.append(L2.dot(ListReducedOrderBasis[n]))
        
    ReducedBasis[n,:]=ListReducedOrderBasis[-1] 
end_time = time.time()
execution_time = end_time - start_time
print(f"Temps d'exécution : {execution_time:.6f} secondes")

Temps d'exécution : 0.007447 secondes

print("-----------------------------------")
print("  Reduced basis accuracy           ")
print("-----------------------------------")
### Offline Errors (H^1_0 and L2 norm)

L2compressionErrors=[]
H10compressionErrors=[]
H10=vector_laplace.assemble(FineBasis['u'],nu=1)
for snap in FineSnapshots:
    ExactSolution =snap
    CompressedSolutionU= ExactSolution@(L2@ReducedBasis.transpose()) #snapshot compression
    ReconstructedCompressedSolution = np.dot(CompressedSolutionU, ReducedBasis)
    
    norml2ExactSolution=np.sqrt(ExactSolution@(L2@ExactSolution))
    normh1ExactSolution=np.sqrt(ExactSolution@(H10@ExactSolution))
    t=np.abs(ReconstructedCompressedSolution-ExactSolution)
    
    if norml2ExactSolution !=0 and normh1ExactSolution != 0:
        relL2Error=np.sqrt(t@L2@t)/norml2ExactSolution
        relH10Error=np.sqrt(t@H10@t)/normh1ExactSolution
    else:
        relL2Error = np.sqrt(t@L2@t) #np.linalg.norm(ReconstructedCompressedSolution-ExactSolution)
        relH10Error=np.sqrt(t@H10@t)
    L2compressionErrors.append(relL2Error)
    H10compressionErrors.append(relH10Error)

print("L2 compression error =", L2compressionErrors)
print("H1 compression error =", H10compressionErrors)

-----------------------------------
  Reduced basis accuracy           
-----------------------------------
L2 compression error = [5.529786242990977e-15, 2.0540757465097742e-14, 0.004569484380782256, 9.279618605250094e-15, 1.9543010683101427e-14, 2.121598047069365e-14, 0.0017752632415230997, 0.0020589462893975245, 0.0012227664193758903, 1.5092299465667746e-14]
H1 compression error = [7.466170361331867e-15, 2.612819937191255e-14, 0.006316823621511566, 1.3728021953919329e-14, 2.8581942857475144e-14, 2.9972234706180816e-14, 0.0024665199322780984, 0.0027363713296409435, 0.001581469755596689, 1.9234340294589924e-14]

cost

If the matrix L2 is dense, then the total cost is in $O(N \ \mathcal{N}_h^2\ Ntrain)$ flops ($O(N \ \mathcal{N}_h \ Ntrain)$ for a sparse matrix), where $\mathcal{N}_h$ is the number of degrees of freedom of one snapshot, $Ntrain$ is the number of snapshots and $N$ the number of modes.

Indeed, MatrixBasisProduct is computed in $O(\mathcal{N}_h^2)$ flops for a dense matrix. We then have two imbricated forloops (in $O(N \ Ntrain)$): The computation of $w$ is in $O(\mathcal{N}_h N)$ and its norm is in $O(\mathcal{N}_h^2)$ flops.

This cost is significantly lowered using a posteriori error estimates.

3. rPOD (random POD)

see details in the last section below

""" rSVD on the correlation matrix """
start_time = time.time()
NumberOfModes=6 #tol
print("number of modes: ",NumberOfModes)

# C=U Sigma VT
# see detailed algorithm steps below (section 4)
U, Sigma, VT = randomized_svd(CorrelationMatrix, 
                              n_oversamples=10, #5 or 10 is better
                              n_components=NumberOfModes, 
                              n_iter=5,  # power iteration
                              #power_iteration_normalizer="none", # or QR
                              random_state=42)  # Répétabilité

EigenValues=TotEigenValues
EigenVectors=VT
print("Eigenvalues: ", EigenValues)
#print(EigenVectors.shape)
RIClist=[]
for i in range(len(EigenValues)):
    RIClist.append(1-np.sum(EigenValues[:i])/np.sum(Sigma)) #must be close to 0
    print("RIC: (",i,")", RIClist[i])

## same steps as in Snapshot-POD
ChangeOfBasisMatrix = np.zeros((NumberOfModes,NumberOfSnapshots))
for j in range(NumberOfModes):
    ChangeOfBasisMatrix[j,:] = EigenVectors[j,:]/np.sqrt(EigenValues[j])

ReducedBasis = np.dot(ChangeOfBasisMatrix,FineSnapshots)
#visualize_velocity(mesh,ReducedBasis[2])

end_time = time.time()
execution_time = end_time - start_time
print(f"Temps d'exécution : {execution_time:.6f} secondes")

number of modes:  6
Eigenvalues:  [ 2.26782714e-02  2.92187702e-03  2.23860208e-04  1.86248477e-05
  2.17765951e-06  1.26997357e-07  4.93353135e-09  1.40613247e-10
  1.09557252e-12 -9.61608314e-19]
RIC: ( 0 ) 1.0
RIC: ( 1 ) 0.12252560668312185
RIC: ( 2 ) 0.009471476038941606
RIC: ( 3 ) 0.0008098105883218576
RIC: ( 4 ) 8.917246601403139e-05
RIC: ( 5 ) 4.913819319996016e-06
RIC: ( 6 ) -4.440892098500626e-16
RIC: ( 7 ) -1.9088965608204944e-07
RIC: ( 8 ) -1.9633030556853726e-07
RIC: ( 9 ) -1.9637269565997428e-07
Temps d'exécution : 0.003512 secondes

Now we can check for our reduced basis accuracy.

print("-----------------------------------")
print("  Reduced basis accuracy           ")
print("-----------------------------------")
### Offline Errors (H^1_0 and L2 norm)

L2compressionErrors=[]
H10compressionErrors=[]
H10=vector_laplace.assemble(FineBasis['u'],nu=1)
for snap in FineSnapshots:
    ExactSolution =snap
    CompressedSolutionU= ExactSolution@(L2@ReducedBasis.transpose()) #snapshot compression
    ReconstructedCompressedSolution = np.dot(CompressedSolutionU, ReducedBasis)
    
    norml2ExactSolution=np.sqrt(ExactSolution@(L2@ExactSolution))
    normh1ExactSolution=np.sqrt(ExactSolution@(H10@ExactSolution))
    t=np.abs(ReconstructedCompressedSolution-ExactSolution)
    
    if norml2ExactSolution !=0 and normh1ExactSolution != 0:
        relL2Error=np.sqrt(t@L2@t)/norml2ExactSolution
        relH10Error=np.sqrt(t@H10@t)/normh1ExactSolution
    else:
        relL2Error = np.sqrt(t@L2@t) #np.linalg.norm(ReconstructedCompressedSolution-ExactSolution)
        relH10Error=np.sqrt(t@H10@t)
    L2compressionErrors.append(relL2Error)
    H10compressionErrors.append(relH10Error)

print("L2 compression error =", L2compressionErrors)
print("H1 compression error =", H10compressionErrors)

-----------------------------------
  Reduced basis accuracy           
-----------------------------------
L2 compression error = [1.3815241858210751e-05, 0.00010641993810958985, 0.00043833388317143893, 0.001211154256116164, 0.0015620109256884352, 0.00040263671234692336, 0.0008616253925896097, 0.00015171542313003271, 0.0006082986187325006, 0.00024666468846097097]
H1 compression error = [1.9399567725164357e-05, 0.0001499851246557588, 0.0006436539912023722, 0.0018634474193020337, 0.0024692938629250977, 0.0006456797464365838, 0.0012552738118668722, 0.00023166001503718088, 0.000837517998819989, 0.00033131133634022204]

cost

If the matrix L2 is dense, then the cost of the correlation matrix computation is in $O(\mathcal{N}_h^2\ Ntrain^2)$ flops where $\mathcal{N}_h$ is the number of degrees of freedom of one snapshot and $Ntrain$ is the number of snapshots.

Then, the dominant cost in the random SVD on the correlation matrix $C$ is in $O (Ntrain^2 (N+p))$ flops (cost of the multiplication of $C \in \mathbb{R}^{Ntrain \times Ntrain}$ by a gaussian matrix $\Omega \in \mathbb{R}^{Ntrain \times (N+p)}$).

4. Truncated SVD

SVD and POD are very close algorithms:

The singular value decomposition (SVD) is widely used in low-rank approximation. It consists in factorizing a matrix $\mathbf{A} \in \mathbb{R}^{m \times n}$ under the form $\mathbf{A} = \mathbf{U} \Sigma \mathbf{V}^T$ where $\mathbf{U} \in \mathbb{R}^{m \times m}$ and $\mathbf{V} \in \mathbb{R}^{n \times n}$ are unitary matrices with $\Sigma = diag(\sigma_1,...,\sigma_p) \in \mathbb{R}^{m \times n}$, and $ \sigma_1 \geq \sigma_2 \geq ··· \geq \sigma_p \geq 0$, with $p = $min$(m,n)$.

The $(\sigma_i)_{i=1,...,p}$ are called singular values of $\mathbf{A}$, $\mathbf{U}=[\Phi_1,\dots,\Phi_m]$ is made up of the left singular vectors of $\mathbf{A}$, whereas the elements of $\mathbf{V}=[\Psi_1,\dots,\Psi_m]$ are its right singular vectors, such that: $\mathbf{A} \Phi_i=\sigma \Psi_i$.

The rank of a diagonal matrix is equal to the number of non zero diagonal terms, and rank($\mathbf{A}$)=rank($\Sigma$), such that if $\mathbf{A}$ has $r$ positive singular values, then rank($\mathbf{A}$) = $r$. We deduce that $\mathbf{A}$ can be written as the following sum:

$\mathbf{A} = \sum_{i=1}^r \sigma_i \Phi_i \Psi_i^T$.

In the POD, $\mathbf{A}^T \mathbf{A}$ is the correlation matrix and we use the fact that $\mathbf{A}^T \mathbf{A}= \mathbf{V} \Sigma^T \Sigma \mathbf{V}^T$.

## SVD on the snapshots 
start_time = time.time()
print(np.array(FineSnapshots).transpose().shape)
U, Sigma, VT = np.linalg.svd(np.array(FineSnapshots).transpose(), full_matrices=False)

# Select the first modes
ReducedBasis = U[:, :NumberOfModes]  # R^(Ndof x NumberOfModes)

# Eigenvalues
EigenValues = Sigma[:NumberOfModes] ** 2  

print("SVD Orthogonality:")
print(ReducedBasis.T @ ReducedBasis)

# RIC
RIC = 1 - np.sum(EigenValues) / np.sum(Sigma ** 2)
print("Relative Information Content (must be close to 0):", RIC)
end_time = time.time()
execution_time = end_time - start_time
print(f"Temps d'exécution : {execution_time:.6f} secondes")

(2178, 10)
SVD Orthogonality:
[[ 1.00000000e+00  1.50920942e-16  7.45931095e-17  1.56125113e-17
  -9.97465999e-18  4.16333634e-17]
 [ 1.50920942e-16  1.00000000e+00  1.25767452e-16  5.20417043e-17
   3.81639165e-17  2.42861287e-17]
 [ 7.45931095e-17  1.25767452e-16  1.00000000e+00  3.66026653e-16
   3.81639165e-17 -3.88686479e-17]
 [ 1.56125113e-17  5.20417043e-17  3.66026653e-16  1.00000000e+00
  -2.18575158e-16 -4.85722573e-17]
 [-9.97465999e-18  3.81639165e-17  3.81639165e-17 -2.18575158e-16
   1.00000000e+00  3.20923843e-17]
 [ 4.16333634e-17  2.42861287e-17 -3.88686479e-17 -4.85722573e-17
   3.20923843e-17  1.00000000e+00]]
Relative Information Content (must be close to 0): 2.0033294889643116e-07
Temps d'exécution : 0.003572 secondes

print("-----------------------------------")
print("  Reduced basis accuracy           ")
print("-----------------------------------")
### Offline Errors (H^1_0 and L2 norm)

L2compressionErrors=[]
H10compressionErrors=[]
H10=vector_laplace.assemble(FineBasis['u'],nu=1)
for snap in FineSnapshots:
    ExactSolution =snap
    CompressedSolutionU= ReducedBasis.transpose()@ExactSolution #snapshot compression
    ReconstructedCompressedSolution = np.dot(CompressedSolutionU, ReducedBasis.transpose())
    
    norml2ExactSolution=np.sqrt(ExactSolution@(L2@ExactSolution))
    normh1ExactSolution=np.sqrt(ExactSolution@(H10@ExactSolution))
    t=np.abs(ReconstructedCompressedSolution-ExactSolution)
    
    if norml2ExactSolution !=0 and normh1ExactSolution != 0:
        relL2Error=np.sqrt(t@L2@t)/norml2ExactSolution
        relH10Error=np.sqrt(t@H10@t)/normh1ExactSolution
    else:
        relL2Error = np.sqrt(t@L2@t) #np.linalg.norm(ReconstructedCompressedSolution-ExactSolution)
        relH10Error=np.sqrt(t@H10@t)
    L2compressionErrors.append(relL2Error)
    H10compressionErrors.append(relH10Error)

print("L2 compression error =", L2compressionErrors)
print("H1 compression error =", H10compressionErrors)

-----------------------------------
  Reduced basis accuracy           
-----------------------------------
L2 compression error = [1.3839653232122614e-05, 0.00010655139750967699, 0.00043876659941552253, 0.0012119568854371354, 0.0015619237548005866, 0.0004033969722890684, 0.0008613698384691759, 0.00015201310598890686, 0.0006078993200356156, 0.0002463245710069301]
H1 compression error = [1.941939702097461e-05, 0.00015005819467710422, 0.0006437531425080919, 0.0018633145654439185, 0.002466401029440625, 0.0006467412928887584, 0.001254080490147989, 0.00023202553425037264, 0.0008359907910302726, 0.0003306576959726763]

cost

The dominant cost in the SVD on the snapshot matrix $A$ is in $O (N_h \times Ntrain^2)$ flops.

5. Random SVD (rSVD)

## rSVD on the snapshots
start_time=time.time()
U, Sigma, VT = randomized_svd(np.array(FineSnapshots).transpose(), 
                              n_components=NumberOfModes, 
                              n_iter=5, 
                              random_state=42) 

ReducedBasis = U[:, :NumberOfModes]  # (Ndof x NumberOfModes)
EigenValues = Sigma[:NumberOfModes] ** 2  

print("rSVD Orthogonality:")
print(ReducedBasis.T @ ReducedBasis)  # Devrait être ≈ Identité
end_time = time.time()
execution_time = end_time - start_time
print(f"Temps d'exécution : {execution_time:.6f} secondes")

rSVD Orthogonality:
[[ 1.00000000e+00 -1.21430643e-17 -9.02056208e-17 -3.29597460e-17
  -3.51281504e-17 -3.81639165e-17]
 [-1.21430643e-17  1.00000000e+00 -9.02056208e-17  1.04083409e-17
  -5.89805982e-17  6.59194921e-17]
 [-9.02056208e-17 -9.02056208e-17  1.00000000e+00  1.26634814e-16
   6.93889390e-18  1.75098651e-17]
 [-3.29597460e-17  1.04083409e-17  1.26634814e-16  1.00000000e+00
   0.00000000e+00 -6.67868538e-17]
 [-3.51281504e-17 -5.89805982e-17  6.93889390e-18  0.00000000e+00
   1.00000000e+00  6.15826834e-17]
 [-3.81639165e-17  6.59194921e-17  1.75098651e-17 -6.67868538e-17
   6.15826834e-17  1.00000000e+00]]
Temps d'exécution : 0.009150 secondes

print("-----------------------------------")
print("  Reduced basis accuracy           ")
print("-----------------------------------")
### Offline Errors (H^1_0 and L2 norm)

L2compressionErrors=[]
H10compressionErrors=[]
H10=vector_laplace.assemble(FineBasis['u'],nu=1)
for snap in FineSnapshots:
    ExactSolution =snap
    CompressedSolutionU= ReducedBasis.transpose()@ExactSolution #snapshot compression
    ReconstructedCompressedSolution = np.dot(CompressedSolutionU, ReducedBasis.transpose())
    
    norml2ExactSolution=np.sqrt(ExactSolution@(L2@ExactSolution))
    normh1ExactSolution=np.sqrt(ExactSolution@(H10@ExactSolution))
    t=np.abs(ReconstructedCompressedSolution-ExactSolution)
    
    if norml2ExactSolution !=0 and normh1ExactSolution != 0:
        relL2Error=np.sqrt(t@L2@t)/norml2ExactSolution
        relH10Error=np.sqrt(t@H10@t)/normh1ExactSolution
    else:
        relL2Error = np.sqrt(t@L2@t) #np.linalg.norm(ReconstructedCompressedSolution-ExactSolution)
        relH10Error=np.sqrt(t@H10@t)
    L2compressionErrors.append(relL2Error)
    H10compressionErrors.append(relH10Error)

print("L2 compression error =", L2compressionErrors)
print("H1 compression error =", H10compressionErrors)

-----------------------------------
  Reduced basis accuracy           
-----------------------------------
L2 compression error = [1.3839653232129974e-05, 0.00010655139750966657, 0.00043876659941556253, 0.001211956885437101, 0.0015619237548005456, 0.0004033969722890735, 0.0008613698384691697, 0.00015201310598888105, 0.0006078993200356165, 0.0002463245710069679]
H1 compression error = [1.941939702088054e-05, 0.00015005819467713285, 0.00064375314250792, 0.001863314565444008, 0.00246640102944091, 0.0006467412928888026, 0.0012540804901480403, 0.00023202553425050946, 0.0008359907910302515, 0.00033065769597253166]

## DETAILS on rSVD:
start_time=time.time()
# if q is small (1 or 2) no power_iteration_normalizer is required but otherwise add QR decomposition at each step (or LU...)
## rSVD on the snapshots

S=np.array(FineSnapshots).transpose()

rng = np.random.default_rng(42)  # for reproductibility purposes
Omega = rng.standard_normal(size=(NumberOfSnapshots, NumberOfModes))  # random matrix #coud also add oversamples

# First we project the snapshots on a smaller random space of size Nt*N
# Step 1 : random projection
#Y = (S@S.T)@(S@S.T)@S@ Omega  

# or with Power iterations to enhance the precision of the approximation of dominant singular vectors (with orthonormalization at each step).
q=5
Y=S@Omega #init
for i in range(q):
    Y = S @ (S.T @ Y)  # Itérations de puissance
    if q>2: 
        Q, _ = np.linalg.qr(Y)  # QR Orthonormalization 
        Y = Q 

# Step 2 : Orthonormalization (QR)
Q, _ = np.linalg.qr(Y)  # (m, k) #S=Q Q^T S

# Step 3  : Reduced projection
B = Q.T @ S  # (k, N)

# Step 4 : SVD on B
U_tilde, S, VT = np.linalg.svd(B, full_matrices=False)
ReducedBasis=Q@U_tilde

end_time = time.time()
execution_time = end_time - start_time
print(f"Temps d'exécution : {execution_time:.6f} secondes")

Temps d'exécution : 0.003637 secondes

print("-----------------------------------")
print("  Reduced basis accuracy           ")
print("-----------------------------------")
### Offline Errors (H^1_0 and L2 norm)

L2compressionErrors=[]
H10compressionErrors=[]
H10=vector_laplace.assemble(FineBasis['u'],nu=1)
for snap in FineSnapshots:
    ExactSolution =snap
    CompressedSolutionU= ReducedBasis.transpose()@ExactSolution #snapshot compression
    ReconstructedCompressedSolution = np.dot(CompressedSolutionU, ReducedBasis.transpose())
    
    norml2ExactSolution=np.sqrt(ExactSolution@(L2@ExactSolution))
    normh1ExactSolution=np.sqrt(ExactSolution@(H10@ExactSolution))
    t=np.abs(ReconstructedCompressedSolution-ExactSolution)
    
    if norml2ExactSolution !=0 and normh1ExactSolution != 0:
        relL2Error=np.sqrt(t@L2@t)/norml2ExactSolution
        relH10Error=np.sqrt(t@H10@t)/normh1ExactSolution
    else:
        relL2Error = np.sqrt(t@L2@t) #np.linalg.norm(ReconstructedCompressedSolution-ExactSolution)
        relH10Error=np.sqrt(t@H10@t)
    L2compressionErrors.append(relL2Error)
    H10compressionErrors.append(relH10Error)

print("L2 compression error =", L2compressionErrors)
print("H1 compression error =", H10compressionErrors)

-----------------------------------
  Reduced basis accuracy           
-----------------------------------
L2 compression error = [1.3839650641800041e-05, 0.00010655138025831539, 0.00043876654776895215, 0.0012119568161507801, 0.0015619237983524252, 0.00040339690992188625, 0.0008613698745842613, 0.00015201308612810389, 0.0006078993568913044, 0.0002463245965823324]
H1 compression error = [1.9419393408006165e-05, 0.0001500581703820479, 0.0006437530670056677, 0.0018633144607903833, 0.002466401103775472, 0.0006467412024139285, 0.0012540805465273642, 0.00023202550603323933, 0.0008359908438539881, 0.0003306577308801652]

cost

The dominant cost in the random SVD on the snapshot matrix $A$ is in $O (N \times Ntrain^2)$ flops.

Authors