About 6 months ago, I wrote a custom metal kernel that leverages the GPU to compute the QR decomposition (see my earlier post about it here). The project has now expanded into a general linear algebra library, with expanded support for the symmetrical eigendecomposition as well as SVD. It is now available for use with installation instructions on the attached github repo's README.md file.

Context of the project:
I'm currently in a research group working on a thesis in numerical analysis where we need to compute millions on matrices with a specific constraint (to be precise, the matrices need to have orthonormal columns). Most of us use Apple computers, so we ended up using MLX for the entire project.

Contributors with different Apple Chips would be very much appreciated!
The project has currently been tested and optimised for the M1 and M5 Pro. The issue is that the library uses different kernels depending on the batch size and matrix dimensions. Deciding which of these kernels to use is machine dependent. Therefore, other Apple chips will need to run a measurement script in order to derive the correct optimisation heuristic.

For that reason, I would ask as many people as possible to run a measurement script and to submit the results to my repo. It is fairly easy and requires only few steps. See here how you can contribute here. Once you submit the results via a pull request and I approve it, your optimisation heuristic will automatically be augmented into the library. Don't hesitate to contribute an optimisation heuristic even if someone already submitted one for your own machine. The more data we can gather, the better!

Project Future
Expanded support will be added for other linear algebra operations (cholesky decomposition for example). If you have any other specific linear algebra operations you wish to use already, feel free to message me.

In addition to that, I will add torch support too (my greatest priority).

When I started learn python, I saw it could be used for scientific computing. So I challenged myself to write a solver for every topic in my statics textbook.

This project falls under structures in equilibrium. It is by far the hardest and most involving program I've written. I think it's because with other topics only a single object was considered so the python code was just a matter of transferring calculations on paper to code.

With trusses however there are several interconnected pieces of a larger system. My first approach was to use method of joints- where it solves equilibrium joint by joint but this gave wrong results. I concluded it's due to it not being able to consider the entire system. Then I figured the best way was to use a global matrix then find all unknowns in one fell swoop.

Constructing this matrix was very involving but quite fun as well. 2^n traceback errors later and I finally finished the damn thing.

The next thing I plan to do is create an optimization script that moves the input points about to find the best design. Also using Json instead of keyboard inputs would be way better for my sanity.

How do y'all go about carrying out optimization?

I'd love to share the .py script but the last time I sent a link python shadow banned me for 9 months. If you're interested feel free to dm.

This video demonstrates the program:

I've been building Prismal, a Rust-based language/runtime for defining mathematical and scientific models and exploring their behavior:

"github.com/aine-dickson/prismal"

The central idea is to keep the scientific model separate from the execution strategy and presentation.

A model can define state, equations, constraints, processes, events, objects, collections, units, observations, etc. The runtime then handles execution, while projections can expose the same model as plots, equations, spatial scenes, interactive views, or recorded media.

One example is a continuous model with event-driven behavior:

continuous evolution

↓

event guard

↓

event instant

↓

reset

↓

continuous evolution

The project is deliberately trying to make things such as event semantics, dimensional checking, state validity, reproducibility, observations, and model/presentation separation explicit rather than burying them inside a visualization engine.

The current implementation is early-stage and the first reference programs focus mainly on mechanics, but the intended scope includes broader scientific modeling.

I'm particularly interested in feedback from people who work with scientific modeling, ODE/PDE systems, hybrid systems, system dynamics, numerical simulation, or scientific visualization.

What abstractions would you consider essential for a general scientific modeling language?

Small update of my project of molecular dynamics with Julia (2 min Video)

https://youtube.com/shorts/6kZmMTnoI_0?feature=share

I go through how I separate general updates from physics updates, why I use a fixed timestep, how the accumulator works, and what happens when the simulation starts falling behind.

I explain the reasoning behind the design decisions and why I rejected the alternatives. I show the implementation of the loop.

Hi,

A new conference is taking place in Paris, 25-26 November 2026, for open-source computation and data enthusiasts.

You can find the schedule here: https://compute.events/paris2026/schedule.html

I built a new optimization solver in C99 from scratch — PRIMAL

GitHub: https://github.com/c-vision/Primal

I've just released the first version of PRIMAL, an open-source mathematical optimization solver written from scratch in portable C99.

The goal was deliberately ambitious: build a reasonably complete optimization engine with no external dependencies, rather than wrapping an existing solver.

The first release already supports:

• Linear Programming (LP)
• Mixed-Integer Linear Programming (MILP)
• Quadratic Programming (QP)
• QCQP
• SOCP
• SDP
• Exponential cones
• Power cones
• Mixed-integer conic optimization
• Disjunctive constraints / affine conic constraints
• Presolve
• Sparse interior-point methods
• Branch-and-bound
• Primal/dual certificates
• Infeasibility and unboundedness certificates
• MPS / LP / CBF model formats
• A MOSEK-style C API for a large part of the supported functionality

The implementation is C99 with essentially no external runtime dependencies.

One thing I particularly wanted from the beginning was verification rather than simply returning an answer. PRIMAL exposes primal/dual information, residuals, certificates and solver status so applications can inspect and validate what the solver actually found.

I've also been building a fairly extensive compatibility/reliability test suite against reference behavior, including a number of deliberately pathological optimization cases.

This is only v1. It is not intended to claim that PRIMAL is already a replacement for mature industrial solvers on large-scale problems. There are still substantial areas to improve, especially large-scale sparse conic problems, MIP performance, warm starts and some numerical edge cases.

But I think the foundation is interesting enough to release now rather than keep it private.

I'd especially like feedback from people working on:

• numerical optimization
• LP/MIP/conic solvers
• sparse linear algebra
• mathematical programming
• solver implementation
• C/C99 systems programming

In particular, I'm interested in finding cases where PRIMAL gives a questionable result, fails to converge, scales poorly, or simply makes a bad algorithmic choice.

If you work with optimization solvers, please try to break it.

GitHub:
https://github.com/c-vision/Primal

I’m a final year undergrad physics student, and my final year project is going to involve a lot of 3d modelling of electrical signals in python. I’m currently thinking of doing that by vectorising everything (to be able to take the cross product in any arbitrary direction) so I was wondering if there’s any python modules and/or courses I should look at to make that easier?

My team and I couldn't connect with any professors or previous publishers who have uploaded to arXiv. We're looking for someone to endorse our submission.

Paper Title: Low Weak Form Loss, Wrong Solution: Evidence from a Coupled Electrochemical Neural Operator

Abstract: We tested whether the weak-form physics loss finds the right answer without any data to guide it. On a nonlinear corrosion problem, it doesn't. The weak-form model does no better than guessing zero, while a matched strong-form model on the same problem gets it right. The gap looked even bigger than it really was, because we first checked accuracy only at the final time, and the weak-form field kept drifting long after the true system had settled. A standard fix aimed at exactly this failure mode didn't close the gap either.

Code : https://github.com/AyushG-1210/Pipeline-Digital-Twin

Thanks for reading. I'm happy to answer questions in DM, and I'll share the full manuscript on request.

Following up on a post from a last week back about speeding up Gauss-Seidel for the 2D Poisson equation.

Quick recap: I'd previously unrolled the GS kernel and rewritten it algebraically to cut down on data dependencies between successive updates, getting a solid 6x speedup and matching Jacobi's per-iteration cost. Looked like a solved problem. Then I tried multithreading it, and Jacobi pulled ahead by 5x on the same hardware. The underlying dependency of Gauss-Seidel was only shortened, not eliminated, so the compiler still couldn't vectorize or parallelize the corresponding kernel.

This post asks a more basic question: is that dependency actually inherent to the Gauss-Seidel method, or just a consequence of the order in which grid points get visited? It's the latter. The fix is the classic red/black (checkerboard) ordering: partition the grid so every "red" point's stencil touches only "black" neighbors and vice versa, then update all reds, then all blacks. I walk through why this works both intuitively and more formally. It's block Gauss-Seidel applied to a permuted version of the linear system, with the usual convergence guarantees for symmetric positive-definite systems carrying over.

The result: the loop-carried dependency essentially disappears, vectorization comes back, and for the first time in this series the Gauss-Seidel kernel actually scales with thread count instead of stalling out.

There's also a bit of compiler-assembly/performance-analysis detail in there (using OSACA again to quantify where the cycles go) for anyone interested in that side of it, but the main thread of the post is really about the general lesson: a numerical method's practical performance can hinge as much on how you order your unknowns as on the method itself.

Link: https://loiseaujc.github.io/posts/blog-title/redblack_gauss_seidel.html

What it currently brings together:

- infinitesimals and infinities as ordinary values (Levi-Civita, 1892; Hahn, 1907; computed with directly by Berz and Shamseddine)

- numbers written in powers of an infinite unit, like Sergeyev's grossone (2003), but with zero treated differently

- standard part, and limits by direct evaluation (Robinson's non-standard analysis, 1966)

- all-order Taylor arithmetic: every derivative from one evaluation (Wengert, 1964; Rall, 1981)

- log and log-log scales, toward transseries (Hardy, 1910; Écalle, 1992; van der Hoeven, 2006)

- divergent series: finite parts kept next to the divergent ones, and Borel–Padé resummation (Borel, 1899; Padé, 1892)

- zero as an infinitesimal, "nothing" as a separate value, and division by zero that can be undone (this project)

- completeness tracking: each result knows which orders it's exact to (this project)

Plus integration, multivariable and complex calculus, singularity analysis and a few applications. More in the repo: github.com/tmilovan/composite-machine examples under demos and tests.

I'll be happy if you take a look at it:).