Power Method Neural Network

Neural Networks

The Power Method Neural Network (PMNN) combines the power iteration method with neural networks to solve linear eigenvalue problems. This approach uses automatic differentiation to compute differential operators and learns eigenfunctions through neural network optimization.

Formulation

For eigenvalue problems of the form $Lu = λ u$ with boundary conditions, PMNN uses a trial function:

u (x) = N_{θ} (x) \cdot ϕ (x)

where $N_{θ} (x)$ is a neural network parameterized by $θ$ , and $ϕ (x)$ is chosen to satisfy boundary conditions. For 1D problems with zero boundary conditions:

ϕ (x) = x (x - L)

which ensures $ϕ (0) = 0$ and $ϕ (L) = 0$ .

PMNN Algorithm

The PMNN algorithm for finding the largest positive eigenvalue:

Build a training dataset $S$ using random sampling
Initialize a neural network with random parameters
For each epoch $k$ :
- Compute $U_{k - 1} (x_{i}) = N_{θ} (x_{i})$ for all $x_{i} \in S$
- Compute $P_{k} = L U_{k - 1}$ using automatic differentiation
- Normalize: $U_{k} = \frac{P _{k}}{∥ P _{k} ∥}$
- Calculate the loss: $loss_{PMNN} (θ) = \frac{1}{N} i = 1 \sum N [U_{k - 1} (x_{i}) - U_{k} (x_{i})]^{2}$
- Update neural network parameters using an optimizer
- If $loss_{PMNN} < ε$ , compute the eigenvalue: $λ = \frac{⟨ L U _{k - 1} , U _{k - 1} ⟩}{⟨ U _{k - 1} , U _{k - 1} ⟩}$ and exit

Implementation

Here is a Python implementation using PyTorch:

import torch
import torch.nn as nn
import torch.optim as optim
import matplotlib.pyplot as plt
import numpy as np

# Neural network definition
class DeepNet(nn.Module):
    def __init__(self, input_dim, hidden_dim, output_dim, num_layers):
        super(DeepNet, self).__init__()
        layers = []
        layers.append(nn.Linear(input_dim, hidden_dim))
        layers.append(nn.Tanh())
        for _ in range(num_layers - 2):
            layers.append(nn.Linear(hidden_dim, hidden_dim))
            layers.append(nn.Tanh())
        layers.append(nn.Linear(hidden_dim, output_dim))
        self.network = nn.Sequential(*layers)
        self.initialize_weights()

    def initialize_weights(self):
        for layer in self.network:
            if isinstance(layer, nn.Linear):
                nn.init.xavier_uniform_(layer.weight)

    def forward(self, x):
        return self.network(x)

# Hyperparameters
L = 1  # Length of region
input_dim = 1
hidden_dim = 20
output_dim = 1
num_layers = 5
N = 1000  # Number of points
N_epoch = 40000
lr = 0.001 
epsilon = 1e-6

# Training dataset
x = torch.linspace(0, L, N).view(-1, 1).requires_grad_(True)

# Neural network initialization
net = DeepNet(input_dim, hidden_dim, output_dim, num_layers)
optimizer = optim.Adam(net.parameters(), lr=lr)

# Boundary conditioning function
def phi(x):
    return x * (L - x)

# Laplacian operator
def laplacian(u, x):
    u_x = torch.autograd.grad(u.sum(), x, create_graph=True, retain_graph=True)[0]
    u_xx = torch.autograd.grad(u_x.sum(), x, create_graph=True, retain_graph=True)[0]
    return u_xx + 100.0*u

# PMNN Loss function
def pmnn_loss(U_k_minus_1, U_k):
    return ((U_k_minus_1 - U_k) ** 2).mean()

# Training loop
U_k_minus_1 = phi(x) * net(x)  # Initial guess
for epoch in range(N_epoch):
    # Compute Laplacian of U_k_minus_1
    P_k = laplacian(U_k_minus_1, x)

    # Normalize P_k to obtain U_k
    norm_P_k = torch.sqrt((P_k ** 2).mean())
    U_k = P_k / norm_P_k

    # Compute PMNN loss
    loss = pmnn_loss(U_k_minus_1, U_k)

    # Backpropagation and optimization
    optimizer.zero_grad()
    loss.backward()
    optimizer.step()

    # Update U_k_minus_1 for the next iteration
    U_k_minus_1 = phi(x) * net(x)

    # Print progress
    if epoch % 1000 == 0:
        print("Epoch {}, Loss: {:.6f}".format(epoch, loss.item()))

    # Stopping criterion
    if loss.item() < epsilon:
        print("Converged at epoch {}, Loss: {:.6f}".format(epoch, loss.item()))
        break

# Compute eigenvalue
x_np = x.detach().numpy().flatten()
u_nn = U_k.detach().numpy().flatten()
LU_k_minus_1_nn = P_k.detach().numpy().flatten()
U_k_minus_1_nn = U_k_minus_1.detach().numpy().flatten()
eigenval = np.dot(LU_k_minus_1_nn, U_k_minus_1_nn)/np.dot(U_k_minus_1_nn, U_k_minus_1_nn)
print("Eigenvalue: ", eigenval)

Key Features

Automatic Differentiation: Uses PyTorch's automatic differentiation to compute differential operators
Boundary Conditions: Incorporated through the trial function $ϕ (x)$
Normalization: Maintains unit norm throughout iteration
Convergence: Stops when the loss falls below a threshold

Applications

Finding dominant eigenvalues of differential operators
Solving eigenvalue problems in quantum mechanics
Computational physics and engineering
High-dimensional eigenvalue problems where traditional methods are expensive

Advantages

Can handle high-dimensional problems efficiently
No need to discretize the domain explicitly
Automatic handling of differential operators
Flexible boundary condition incorporation

References

Yang, Q., Deng, Y., Yang, Y., He, Q., & Zhang, S. (2023). Neural networks based on power method and inverse power method for solving linear eigenvalue problems. Computers & Mathematics with Applications, 147, 14-24.