[RFC]: add linear algebra functionality

[dhruvan2006]

Full name

Dhruvan Gnanadhandayuthapani

University status

Yes

University name

Delft University of Technology

University program

B.Sc. Computer Science and Engineering

Expected graduation

2026

Short biography

I am Dhruvan and am currently 3rd (final) year student at the Delft University of Technology (TU Delft) in the Netherlands. Professionally, I was an intern at Google Amsterdam during the summer of 2025 where I had to heavily utilize TypeScript to build a full stack application that integrated AI workflows. I have also served as a Teaching Assistant in my university, notably for the course Machine Learning. Outside of this, I have a keen interest for high frequency trading (C++) and general machine learning (Python/NumPy) which are evidenced by my projects below.

Timezone

UTC+01:00

Contact details

email:dhruvan2006@gmail.com,github:dhruvan2006

Platform

Windows

Editor

I primarily use Neovim on WSL (for a UNIX environment). I find it very lightweight and customizable and am able to find things and move between files very quickly. It helps me reduce my reliance on AI tools (I do use inline completion) and helps me to reduce distraction and focus on the code.

Programming experience

I have previous experience during my internship at Google building full stack applications using TypeScript. The tech stack composed of Next.js and integrating AI workflows into it. Unfortunately, I am not able to share the code or work with you as it is an internal tool.

My first big project that utilized JavaScript is the following project. I wanted a website that showcases my in-house algorithms and strategies and decided to build it with Next.js. The website is live here and the repo is available here.

I also have a very keen interest for high frequency trading and doing system level optimizations. I primarily use C/C++ for this task and have built a project called OrderBook where I focused on building an ultra-low latency matching engine (sub-100ns median latency).

JavaScript experience

I have primarily used JavaScript through frameworks. During my internship at Google last summer, we extensively used Next.js to craft our full-stack application.

Although Javascript is single-threaded, the existence of the event loop makes writing async/await code very natural. Along with the Promise API, working with any kind of IO in JS is by far the most enjoyable for me.

Coming from strongly typed languages my least favorite feature of JavaScript is type coercion. JavaScript tries to be "helpful" by converting types on the fly but it leads to results that make no sense on the surface:

Node.js experience

While I have used Node.js in server-side environments (as noted above in my work with Next.js and also having used Express.js), my contributions to stdlib (adding C native addons) prompted me to dive into Node’s internals. Specifically, I have been studying the Node-API and how it allows developers to create native addons.

C/Fortran experience

C

I have a keen interest in HFT and C++ is the standard for performance critical applications. I have built highly cache and latency aware systems before (most notably an order book with sub-100ns median latency). I have also spent considerable time studying the source code of other big open source projects like TA-Lib and PostgreSQL.

Fortran

I do not have any prior experience with Fortran outside the scope of stdlib. I have spent the past month studying the syntax of Fortran 77 along with reference implementations of BLAS and LAPACK routines.

For a project like stdlib, my main takeaway is the fact that memory layout differences between Fortran (column-major) and JS/C (row-major) should be handled efficiently along with other optimization techniques like cache blocking and loop unrolling.

Interest in stdlib

I love the rigor in stdlib. The maintainers go to great lengths to make sure that packages follow a consistent style guide and have 100% test coverage, benchmarks and extensive API documentation.

My favorite feature of stdlib is the ndarray. Native JS does not guarantee that arrays are stored as contiguous blocks of memory. stdlib fixes this by backing multidimensional arrays with typed arrays. It allows views and slicing as well as zero copy operations for transposing (this allows optimizations when implementing linalg functions like multiplication of a matrix with its transpose).

Another worthy mention is office hours! I really enjoyed asking questions and interacting with the maintainers. This level human contact is not seen in many other projects.

Version control

Yes

Contributions to stdlib

I started off by studying existing JS and C implementation of BLAS routines in stdlib. This led me to find an incorrect function call stdlib-js/stdlib#10707 (merged). I also noticed an improvement that could be made to dgemm's tests stdlib-js/stdlib#10709.

Moving forward, to better estimate the time requirement for implementing BLAS and LAPACK routines, I went ahead and added the C implementation for dgemm stdlib-js/stdlib#10721 and the JS implementation for dlange stdlib-js/stdlib#10661. I also reviewed stdlib-js/stdlib#7016. All three of these are blockers for my Linear Algebra proposal.

Merged

stdlib-js/stdlib#10707

Unmerged

stdlib-js/stdlib#10661
stdlib-js/stdlib#10709
stdlib-js/stdlib#10721

Reviews

stdlib-js/stdlib#7016

stdlib showcase

I wanted to showcase the linear algebra and mathematical capabilities of stdlib. I also wanted to compare the effectiveness of JavaScript as a language for ML. Therefore, I decided to train a Multilayer Perceptron (MLP) on the MNIST dataset. I compared the implementation in JS (stdlib) vs JS (stdlib with native addons) vs Python (Numpy). There is also a live demo of the final model (using UMD modules).

Repo: https://github.com/dhruvan2006/stdlib-mnist/
Live demo: https://dhruvan2006.github.io/stdlib-mnist/

Benchmark results:

Implementation	Total Time (10 Epochs)	Average Time / Epoch	Speedup
Vanilla JS (Baseline)	309.0s	~30.9s	1.00x
NumPy	191.7s	~19.2s	1.61x
Stdlib Addons	173.7s	~17.4s	1.78x

Goals

Abstract

stdlib currently does not have a user friendly API for linear algebra. The library currently only has a subset of BLAS and LAPACK routines implemented. The goal of this proposal is to add a linear algebra namespace to stdlib (@stdlib/linalg/) which provides an intuitive and flexible API. All the underlying BLAS and LAPACK routines will be implemented in both JS and C.

API

Below is the subset of the linear algebra functionality that will be implemented along with a technical description of their logic. This subset has been chosen to be comprehensive enough to serve as a starting point for this library and also such that they have many overlapping dependencies (BLAS and LAPACK)

NOTE: The dependencies for qr are a subset of the ones proposed in #184 and will therefore become a stretch goal if that proposal is not selected.

The routines in red are yet to be implemented. The routines in orange have PRs open and need to be reviewed. The PRs in green have already been implemented.

1. linalg/dot

Computes the dot product or generalized matrix multiplication for scalars, vectors, and matrices.

Spec

• Check if a or b is scalar:
  • If length==1 just multiply
  • If vector, then axpy
  • If matrix, then for loop axpy on longest dimension of matrix
• Both are vectors: Level 1 BLAS *dot
• Matrix vector: Level 2 BLAS gemv
• Vector matrix: Remember that $vA = (A^Tv^T)^T$, therefore Level 2 BLAS gemv
• Matrix matrix: Remember to get correct dimensions
  • If case of matrix times its transpose, use syrk (don't do $O(n^2)$ check, try to only check the flags)
  • Else gemm

2. linalg/matmul

Performs general matrix multiplication for vectors and matrices, with behavior similar to dot but generalized for higher-dimensional arrays.

Spec

• If 1-D or vector, do the same as dot:
  • 1 x n times n x 1: dot
  • Scalar times vec:
  • Vector times matrix: Use gemv with identity $vA = (A^Tv^T)^T$
  • Matrix times vector: Use gemv normally
• Matrix times matrix: Use the identity $A \times B = C \implies B^T \times A^T = C^T$ to use the closer strided column or row
  • In case of matrix times its transpose, use syrk
  • For all other cases, use gemm (gemv cannot conjugate the vector)

3. linalg/norm

Computes vector and matrix norms for a given input.

Spec

ord	matrix	vectors
None	Frobenius norm	2-norm
'fro'	Frobenius norm	-
'nuc'	nuclear norm	-
inf	max(sum(abs(x), axis=1))	max(abs(x))
-inf	min(sum(abs(x), axis=1))	min(abs(x))
0	-	as below
1	max(sum(abs(x), axis=0))	as below
-1	min(sum(abs(x), axis=0))	as below
2	2-norm (largest sing. value)	as below
other	-	sum(abs(x)**ord)**(1./ord)

• Fast path for 1D vector:

  • Frobenius norm or 2-norm: dot product with itself (if complex, sum of real dot and imaginary dot)

• 1D vector norms:

  • If ord == inf, return the absolute maximum value
  • If ord == -inf, return the absolute minimum value
  • If ord == 0, return the count of non-zero elements
  • If ord == 1, return the sum of absolute values
  • If ord == 2, standard Euclidean norm
  • Else generic p norm: $(\sum |x|^p)^{1/p}$

• Matrix norms:

  • (Requires SVD) If ord == 2, use SVD to find the largest singular value
  • (Requires SVD) If ord == -2, use SVD to find the smallest singular value
  • If ord == 1, maximum column sum
  • If ord == -1, minimum column sum
  • If ord == inf, maximum row sum
  • If ord == -inf, minimum row sum
  • If Frobenius: $(\sqrt{\sum |x_{i,j}|^2})$
  • (Requires SVD) If nuclear, use SVD to find the sum of all singular values
  • Else, throw error

4. linalg/cholesky

Computes the Cholesky decomposition of a symmetric (Hermitian) positive-definite matrix.

Spec

Call potrf

5. linalg/lu

Computes the LU decomposition of a matrix.

Spec

Call getrf()

6. linalg/det

Computes the determinant of a square matrix.

Spec

Given LU decomposition:
$A = P^{-1}LU$ of square matrix $A$, the determinant of A
$$det(A) = det(P^{-1})det(L)det(U)$$
$$det(A) = (-1)^S (\prod_{i=1}^{n}l_{ii})(\prod_{i=1}^{n}u_{ii})$$
where $S$ is the number of row exchanges in the decomposition

Remember to handle singular matrices.
Make a call to getrf() to LU factorize. Then use above formula.

7. linalg/inv

Computes the inverse of a square matrix.

Spec

Call gesv() with $AX = B$, where $B$ is the identity matrix $I$. The resulting $X$ is $A^{-1}$.

8. linalg/solve

Solves a system of linear equations of the form: :$AX = B$

Spec

Call gesv()

9. linalg/qr

Computes the QR decomposition of a matrix.

Spec

Let K = min(M, N), then

• 'reduced' : returns Q, R with dimensions (..., M, K), (..., K, N)
• 'complete' : returns Q, R with dimensions (..., M, M), (..., M, N)
• 'r' : returns R only with dimensions (..., K, N)
• 'raw' : returns h, tau with dimensions (..., N, M), (..., K,)

Implementation:

• Run geqrf()
• If mode 'complete', call gqr() for $Q$ with $M \times M$ dimensions
• If mode 'reduced', call gqr() for $Q$ with $M \times min(M, N)$ dimensions

BLAS

The following BLAS routines have to be implemented. Both have open PRs.

• dtrmv (C feat: add c implementation for blas/base/dtrmv stdlib#7016)
• dgemm (C feat: add c implementation for blas/base/dgemm stdlib#10721)

LAPACK

The list of LAPACK routines is more extensive. The routines for QR have lesser priority and overlap with those required for SVD (proposed in #184). They are italicized.

2 JS (+4 for SVD) have PRs open

6 JS (+4 for SVD) routines need to be implemented
10 C (+8 for SVD) routines need to be implemented

• dlarf1f (C)
• dlarfg (JS feat: add lapack/base/dlarfg stdlib#7109, C)
• dtrmm (JS feat: add blas/base/dtrmm stdlib#7366, C)
• dlamch (C)
• dlaswp (C)
• dtrsm (JS feat: add blas/base/dtrsm stdlib#7359, C)
• dgeqr2 (JS, C)
• dlarfb (JS feat: add lapack/base/dlarfb stdlib#7833, C)
• dlarft (JS, C)
• dgetrf2 (JS, C)
• dgetrs (JS, C)
• dpotrf2 (JS, C)
• dlange (JS feat: add lapack/base/dlange stdlib#10661, C)
• dgeqrf (JS, C)
• dgetrf (JS, C)
• dpotrf (JS, C)
• dorgqr (JS, C)
• dorg2r (JS feat: add lapack/base/dorg2r stdlib#7784, C)

Why this project?

I have a keen interest for performance critical computing and machine learning, which drives my desire to bring these tools to the JavaScript ecosystem. My showcase project demonstrates that JS has the capability of matching, and in some cases exceeding, the speed of the the gold standard NumPy.

Despite the sheer scale of the ecosystem (with TypeScript and JavaScript currently ranking as the first and third most popular languages on GitHub), JavaScript still lacks a standardized high performance linear algebra suite. This project would provide the foundational primitives that will help build neural networks, physics engines, robotics simulations, etc. By removing the barrier to entry for millions of web developers, I would like to help democratize AI.

stdlib is the perfect project for this implementation because of the already existing ndarray infrastructure which provides many necessary memory management primitives. The project also has an extensive implementation of BLAS routines and already has a well-defined and seamless interface with C, Fortran and WASM.

Qualifications

I have taken multiple courses throughout my Bachelor in topics related to Linear Algebra. I have completed CSE1205 - Linear Algebra which introduced me to topics such as matrix algebra (sum, product, inverse) solving a linear system of equations, performing LU decomposition, determinants, eigenvalues/vectors and the like.

This theory was put into practice in three other courses that I studied that heavily relied on concepts of Linear Algebra, viz., CSE2510 - Machine Learning, CSE2530 - Computational Intelligence and CSE3130 - Introduction to Quantum Science.

Prior art

The main aim of this library has been implemented very well in other programming languages and has also been attempted before in the JavaScript ecosystem.

SciPy and NumPy (Python)

NumPy and SciPy are the industry standards for numerical computing in Python. They are thin wrappers around high performance C and Fortran libraries. The proposed API for stdlib draws inspiration from numpy.linalg, especially in the definitions of operations for dot, matmul, etc.

Julia

Julia was designed for high performance numerical computing. It achieves near C level performance through JIT compilation while maintaining a high level syntax. Julia’s standard library includes comprehensive linear algebra support backed by optimized BLAS and LAPACK implementations. It serves as a strong example of how a language can provide both usability and performance without relying on a two-language model

scijs

scijs is a modular collection of packages that pioneered multidimensional arrays and numerical computing in JavaScript. While the ecosystem introduced important primitives such as ndarray and low-level linear algebra routines, it is somewhat fragmented and many packages are no longer actively maintained, making it less suitable as a cohesive modern solution.

Math.js

Math.js is a popular and feature rich library for JavaScript. However, it does not prioritize performance and relies heavily on pure JavaScript implementations, making it less suitable for large scale or compute intensive linear algebra tasks.

Tensorflow.js

While the library is primarily for deep learning, tensorflow.js does contain highly optimized linear algebra routines. However, they are mostly targeted towards using the GPU (by utilizing WebGL or WebGPU) and make it less ideal for general purpose CPU bound tasks.

BLAS and LAPACK

Low level numerical computation in most high performance systems is built on standardized libraries such as BLAS (Basic Linear Algebra Subprograms) and LAPACK (Linear Algebra Package). These libraries define efficient and well tested routines for vector and matrix operations. This proposal will primarily use Netlib's reference implementations as the source of truth when implementing BLAS and LAPACK routines.

Commitment

I can consistently devote ~30 hr/week during the entire length of the project, except on the week of June 22 (week 5) where I will have to present my Bachelor thesis. During that specific week, I will remain active but with slightly reduced availability and have planned the Schedule below based on this information.

I will also be devoting time during the community bonding period to review all the open PRs (and take over them if required) mentioned in my proposal which would be blockers during the coding period.

I also noticed many abandoned PRs from previous GSoC participants and will ensure that this does not happen. I will take the time to apply changes from the review and reiterate on the work even after the program ends to make sure that everything that has been promised by me will be delivered.

I am also interested in slowly adding more linear algebra functionality (eigenvalues, least squares, pinv, schur, etc.) outside of those mentioned in the proposal in the following year to achieve the end goal of having a complete linear algebra package in stdlib.

Schedule

The primary goal is to maximize concurrency while maintaining sufficient time for review cycles by leaving a minimum 2 week long buffer between an implementation and its dependency as well as between the JS and C versions of each routine. The secondary goal is to balance the weekly workload by using the line count of each routine as a rough proxy for task “difficulty”.

For each week I will include the BLAS/LAPACK routines and linear algebra packages that I will be working on completing.

This sums up to a total of:

• 8 PR reviews (might need to take some over due to inactivity)
• 6 LAPACK JavaScript routines
• 10 LAPACK C routines
• 9 Linear algebra functions

Assuming a 12 week schedule,

• Community Bonding Period:
  Reviewing (and/or taking over) the following 8 PRs which are blockers to the features I want to implement in the standard coding period:
  feat: add c implementation for blas/base/dtrmv stdlib#7016, feat: add c implementation for blas/base/dgemm stdlib#10721, feat: add lapack/base/dlarfg stdlib#7109, feat: add blas/base/dtrmm stdlib#7366, feat: add blas/base/dtrsm stdlib#7359, feat: add lapack/base/dlarfb stdlib#7833, feat: add lapack/base/dlange stdlib#10661, feat: add lapack/base/dorg2r stdlib#7784).

The numbers in the brackets are the lines of code for each routine to act as a rough metric of the "difficulty" of it.

• Week 1:
  JS: dgetrf2 (155), dgetrs (103)
  C: dlamch (76)

• Week 2:
  JS: dpotrf2 (127)
  C: dlaswp (74), dtrsm (259)

• Week 3:
  JS: dgetrf (117)
  C: dlange (94), dgetrf2 (155)

• Week 4:
  JS: dpotrf (140)
  C: dgetrs (103), dpotrf2 (127)

• Week 5:
  JS: dgesv (53)
  LA: norm (medium)

• Week 6: (midterm)
  Midterm submission!
  C: dgetrf (117), dpotrf (140)
  LA: norm (medium; contd from previous week)

norm is continued from previous week because its the first LA routine and might cause delays to get the boilerplate up and running (API style, new kind of documentation, etc.)

• Week 7:
  C: dgesv (53)
  LA: matmul (hard), dot (hard)

• Week 8:
  LA: solve (hard)

• Week 9:
  LA: lu (medium), det (medium)

• Week 10:
  LA: inv (medium), cholesky (medium)

• Week 11:
  LA: qr (hard) (if dependencies are met, if not can provide a helping hand to implement dependencies)

• Week 12: Catch all week for any delays in the above and/or any unforeseen circumstances
  LA: qr (hard) (if dependencies are met, if not can provide a helping hand to implement dependencies)

• Final Week:

  • Make sure that everything is completed and merged.
  • Project submission!

Notes:

• The community bonding period is a 3 week period built into GSoC to help you get to know the project community and participate in project discussion. This is an opportunity for you to setup your local development environment, learn how the project's source control works, refine your project plan, read any necessary documentation, and otherwise prepare to execute on your project project proposal.
• Usually, even week 1 deliverables include some code.
• By week 6, you need enough done at this point for your mentor to evaluate your progress and pass you. Usually, you want to be a bit more than halfway done.
• By week 11, you may want to "code freeze" and focus on completing any tests and/or documentation.
• During the final week, you'll be submitting your project.

Related issues

The GSoC idea: #28

Checklist

• I have read and understood the Code of Conduct.
• I have read and understood the application materials found in this repository.
• I understand that plagiarism will not be tolerated, and I have authored this application in my own words.
• I have read and understood the patch requirement which is necessary for my application to be considered for acceptance.
• I have read and understood the stdlib showcase requirement which is necessary for my application to be considered for acceptance.
• The issue name begins with [RFC]: and succinctly describes your proposal.
• I understand that, in order to apply to be a GSoC contributor, I must submit my final application to https://summerofcode.withgoogle.com/ before the submission deadline.

[kgryte]

@dhruvan2006 Thank you for opening this draft proposal. A few comments:

1. Thank you for your showcase. The live link did not work for me, but I did browse the code. And seeing the code certainly reinforced the need for more ergonomic linear algebra APIs, as you suggest in this proposal.
2. I worry a bit about the scope of your proposal. I see Cholesky, QR, LU, matrix inversion, solving systems of linear equations, SVD, norms, and more. That is quite ambitious, and I am not sure that your dependency analysis fully captures the scope of the required dependency tree.
3. One of the most difficult aspects of implement BLAS/LAPACK routines is testing them. Especially in LAPACK, routines are "tested" based on their usage. Meaning, a subroutine is not explicitly tested, only the top-level routine, with the subroutine being considered "correct" so long as nothing breaks in the top-level routine. As such, in the past, we've had to build our own test suites. This can be a long process, especially as we implement BLIS-style APIs for row-/column-major generality. You should probably touch on your testing strategy.
4. Some of your proposed top-level routines are off in their "spec". E.g., dot should not be overloaded to perform general matrix multiplication. matmul should be strictly for matrix-matrix multiplication, etc.
5. One of the most glaring omissions for building top-level APIs is how we go from ndarrays to strided interfaces. This is an ongoing story, as seen in the stats/* namespace where we now have the ability to perform reductions which call into strided routines. Similarly in blas/ext. We have yet to build the machinery for calling into >1D strided routines and certainly haven't built machinery for strided routines accepting arrays of mixed rank. There will likely need to be R&D to figure out how to create this bridge between arbitrary views on multi-dimensional data and the underlying strided routines which do the actual computation. Until that has been completed, top-level APIs will be hard to implement.