Узнайте, как построить профессиональную вычислительную платформу на чистом Python без NumPy и SciPy. Реализуйте линейную алгебру, квантовые симуляции и нейросети с нуля. Попробуйте сами!
Вы когда-нибудь задумывались, что скрывается под капотом NumPy, SciPy или PyTorch? Большинство разработчиков воспринимают эти библиотеки как чёрную магию. А что, если попробовать построить профессиональную вычислительную платформу, используя только стандартную библиотеку Python? Это не просто упражнение — это глубокое погружение в алгоритмическую оптимизацию, управление памятью и чистую математику. В этой статье я расскажу, как мне удалось создать такую платформу, и покажу ключевые фрагменты кода.
Наивная реализация определителя матрицы через разложение по строке имеет факториальную сложность — O(N!). Для матрицы 10×10 это уже миллионы операций. Я заменил её на LU-разложение с частичным выбором главного элемента, которое работает за O(N³). Вот как выглядит ключевая часть:
def lu_decomposition(A):
"""LU-разложение с частичным выбором главного элемента.
Возвращает P, L, U такие, что P @ A = L @ U.
"""
n = len(A)
L = [[0.0] * n for _ in range(n)]
U = [[0.0] * n for _ in range(n)]
P = list(range(n)) # матрица перестановок
for k in range(n):
# поиск главного элемента
max_row = max(range(k, n), key=lambda i: abs(A[i][k]))
if A[max_row][k] == 0:
raise ValueError("Матрица вырождена")
# перестановка строк
P[k], P[max_row] = P[max_row], P[k]
A[k], A[max_row] = A[max_row], A[k]
L[k][k] = 1.0
for i in range(k + 1, n):
L[i][k] = A[i][k] / A[k][k]
for j in range(k, n):
A[i][j] -= L[i][k] * A[k][j]
for i in range(n):
for j in range(i, n):
U[i][j] = A[i][j]
return P, L, UЭтот код — основа для решения систем линейных уравнений, вычисления обратных матриц и определителей. Без NumPy пришлось вручную управлять индексами и избегать избыточных копирований.
Симуляция квантовых систем обычно требует перемножения матриц размерности 2ⁿ, что быстро становится неподъёмным. Я реализовал вероятностный алгоритм выборки, который работает за O(1) для каждого измерения. Идея в том, чтобы хранить вектор состояния и применять операторы без полного пересчёта:
import random
import math
class QuantumState:
def __init__(self, num_qubits):
self.n = num_qubits
self.amplitudes = [0.0] * (1 << num_qubits)
self.amplitudes[0] = 1.0 # начальное состояние |0...0⟩
def apply_gate(self, gate, qubits):
"""Применение гейта к указанным кубитам.
gate — матрица 2x2, заданная списком списков.
"""
# здесь реализация тензорного произведения
pass
def measure(self):
"""Измерение с коллапсом состояния за O(1)."""
probs = [abs(a)**2 for a in self.amplitudes]
r = random.random()
cumulative = 0.0
for i, p in enumerate(probs):
cumulative += p
if r < cumulative:
# коллапс в состояние |i⟩
self.amplitudes = [0.0] * len(self.amplitudes)
self.amplitudes[i] = 1.0
return i
return len(probs) - 1Алгоритм выборки с накоплением вероятностей позволяет получить результат без многократного пересчёта матриц.
Полноценный многослойный перцептрон (MLP) с кастомным оптимизатором Adam — это вызов. Пришлось реализовать обратное распространение ошибки вручную, без autograd. Вот фрагмент класса нейросети:
import math
class MLP:
def __init__(self, layers):
self.layers = layers # список кортежей (in_features, out_features)
self.weights = []
self.biases = []
for in_dim, out_dim in layers:
# инициализация Xavier
std = math.sqrt(2.0 / (in_dim + out_dim))
w = [[random.gauss(0, std) for _ in range(out_dim)] for _ in range(in_dim)]
b = [0.0] * out_dim
self.weights.append(w)
self.biases.append(b)
def forward(self, x):
activations = [x]
for w, b in zip(self.weights, self.biases):
z = [sum(xi * w[i][j] for i, xi in enumerate(x)) + b[j] for j in range(len(b))]
x = [1.0 / (1.0 + math.exp(-zi)) for zi in z] # сигмоида
activations.append(x)
return x, activations
def backward(self, y_true, activations):
# градиенты для выходного слоя
grad = [(a - t) * a * (1 - a) for a, t in zip(activations[-1], y_true)]
# ... и так для каждого слоя
passОптимизатор Adam сохраняет моменты первого и второго порядка для каждого параметра, что требует аккуратного управления состоянием.
Быстрое преобразование Фурье (FFT) — классический алгоритм. Моя реализация использует дополнение нулями до степени двойки и рекурсивное разделение:
def fft(x):
n = len(x)
if n <= 1:
return x
# дополнение нулями до степени двойки
m = 1
while m < n:
m <<= 1
x = x + [0.0] * (m - n)
n = m
# рекурсивное БПФ
even = fft(x[0::2])
odd = fft(x[1::2])
factor = [cmath.exp(-2j * cmath.pi * k / n) * odd[k] for k in range(n // 2)]
return [even[k] + factor[k] for k in range(n // 2)] + \
[even[k] - factor[k] for k in range(n // 2)]Использование теоремы о свёртке позволяет вычислять свёртку за O(N log N) вместо O(N²).
Движок сам генерирует гипотезы (например, «средние двух выборок различаются») и проверяет их с помощью ANOVA или t-теста Уэлча. Реализация критерия Уэлча:
import math
def welch_t_test(sample1, sample2):
n1, n2 = len(sample1), len(sample2)
mean1 = sum(sample1) / n1
mean2 = sum(sample2) / n2
var1 = sum((x - mean1)**2 for x in sample1) / (n1 - 1)
var2 = sum((x - mean2)**2 for x in sample2) / (n2 - 1)
t = (mean1 - mean2) / math.sqrt(var1/n1 + var2/n2)
df = (var1/n1 + var2/n2)**2 / ((var1/n1)**2/(n1-1) + (var2/n2)**2/(n2-1))
return t, dfЗначение p вычисляется через аппроксимацию распределения Стьюдента.
Чтобы доказать, что это не игрушка, я обеспечил:
Всё это доступно в открытом репозитории.
Не бойтесь заглядывать под капот. Написание собственной реализации алгоритмов — лучший способ понять их по-настоящему. Попробуйте сами: установите pip install cognitive-discovery-platform, изучите код и, возможно, внесите свой вклад. Звёзды и пул-реквесты приветствуются!
Хочешь закрепить знания на практике?
Решай задачи на Algolit — интерактивная платформа для обучения
Начать бесплатно →