diff --git a/matrix.d.ts b/matrix.d.ts
index 82689f4..4e02b5e 100644
--- a/matrix.d.ts
+++ b/matrix.d.ts
@@ -626,6 +626,16 @@ export abstract class AbstractMatrix {
    */
   gram(): Matrix;
 
+  /**
+   * Returns the matrix product between `this` and its transpose (`this · thisᵀ`),
+   * optionally weighting each column by `scale` (`this · diag(scale) · thisᵀ`).
+   * The result is symmetric, so only its upper triangle is computed and mirrored
+   * and the transpose is never materialized, making it about twice as fast as
+   * `this.mmul(this.transpose())`.
+   * @param scale - Optional per-column factors (length equal to the number of columns).
+   */
+  mmulByTranspose(scale?: ArrayLike<number>): Matrix;
+
   /**
    * Returns the square matrix raised to the given power
    * @param scalar - the non-negative integer power to raise this matrix to
diff --git a/src/__tests__/matrix/utility.test.js b/src/__tests__/matrix/utility.test.js
index f686ced..ff0d867 100644
--- a/src/__tests__/matrix/utility.test.js
+++ b/src/__tests__/matrix/utility.test.js
@@ -332,6 +332,45 @@ describe('utility methods', () => {
     ]);
   });
 
+  it('mmulByTranspose', () => {
+    let matrix = new Matrix([
+      [1, 2, 3],
+      [4, 5, 6],
+    ]);
+    expect(matrix.mmulByTranspose().to2DArray()).toStrictEqual([
+      [14, 32],
+      [32, 77],
+    ]);
+    // equivalent to this.mmul(this.transpose()), but symmetric and without transpose
+    expect(matrix.mmulByTranspose().to2DArray()).toStrictEqual(
+      matrix.mmul(matrix.transpose()).to2DArray(),
+    );
+  });
+
+  it('mmulByTranspose with column scale', () => {
+    let matrix = new Matrix([
+      [1, 2, 3],
+      [4, 5, 6],
+    ]);
+    let scale = [2, 1, 0.5];
+
+    // this · diag(scale) · thisᵀ
+    expect(matrix.mmulByTranspose(scale).to2DArray()).toStrictEqual([
+      [2 * 1 + 1 * 4 + 0.5 * 9, 2 * 4 + 1 * 10 + 0.5 * 18],
+      [2 * 4 + 1 * 10 + 0.5 * 18, 2 * 16 + 1 * 25 + 0.5 * 36],
+    ]);
+
+    // equivalent to this.mmul(diag(scale)).mmul(this.transpose())
+    const diag = Matrix.diag(scale);
+    expect(matrix.mmulByTranspose(scale).to2DArray()).toStrictEqual(
+      matrix.mmul(diag).mmul(matrix.transpose()).to2DArray(),
+    );
+
+    expect(() => matrix.mmulByTranspose([1, 2])).toThrow(
+      'scale must have one value per column',
+    );
+  });
+
   it('mmul strassen on empty matrices', () => {
     // https://github.com/mljs/matrix/issues/114
     // while the mathematically correct result is 0x0, we assert a 2x2 padded result that the current implementation produces
diff --git a/src/matrix.js b/src/matrix.js
index 805df48..c7a5446 100644
--- a/src/matrix.js
+++ b/src/matrix.js
@@ -908,6 +908,44 @@ export class AbstractMatrix {
     return result;
   }
 
+  mmulByTranspose(scale) {
+    let m = this.rows;
+    let n = this.columns;
+
+    if (scale !== undefined && scale.length !== n) {
+      throw new RangeError('scale must have one value per column');
+    }
+
+    let result = new Matrix(m, m);
+
+    // result = this · diag(scale) · thisᵀ is symmetric, so only the upper
+    // triangle is computed and mirrored, and the transpose is never
+    // materialized. `scale` (one factor per column) is folded into one operand.
+    let rowj = new Float64Array(n);
+    for (let j = 0; j < m; j++) {
+      if (scale === undefined) {
+        for (let k = 0; k < n; k++) {
+          rowj[k] = this.get(j, k);
+        }
+      } else {
+        for (let k = 0; k < n; k++) {
+          rowj[k] = scale[k] * this.get(j, k);
+        }
+      }
+
+      for (let i = j; i < m; i++) {
+        let s = 0;
+        for (let k = 0; k < n; k++) {
+          s += this.get(i, k) * rowj[k];
+        }
+
+        result.set(i, j, s);
+        result.set(j, i, s);
+      }
+    }
+    return result;
+  }
+
   mpow(scalar) {
     if (!this.isSquare()) {
       throw new RangeError('Matrix must be square');