ארבעת תתי־המרחבים היסודיים

function transpose(A) {
  return A[0].map((_, j) => A.map((row) => row[j]));
}

function rref(input, relativeTolerance=1e-12) {
  const matrix = input.map((row) => row.map(Number));
  const rowScales = matrix.map((row) => Math.max(...row.map(Math.abs)));
  const pivots = [];
  let pivotRow = 0;
  for (let column = 0;
       column < matrix[0].length && pivotRow < matrix.length;
       column += 1) {
    let best = -1;
    let bestRatio = 0;
    for (let row = pivotRow; row < matrix.length; row += 1) {
      const ratio = rowScales[row] === 0
        ? 0
        : Math.abs(matrix[row][column])/rowScales[row];
      if (ratio > bestRatio) {
        bestRatio = ratio;
        best = row;
      }
    }
    if (best === -1 || bestRatio <= relativeTolerance) continue;
    [matrix[pivotRow], matrix[best]] = [matrix[best], matrix[pivotRow]];
    [rowScales[pivotRow], rowScales[best]] = [rowScales[best], rowScales[pivotRow]];
    const scale = matrix[pivotRow][column];
    matrix[pivotRow] = matrix[pivotRow].map((value) => value/scale);
    for (let row = 0; row < matrix.length; row += 1) {
      if (row === pivotRow) continue;
      const factor = matrix[row][column];
      matrix[row] = matrix[row].map(
        (value, j) => value - factor*matrix[pivotRow][j],
      );
    }
    pivots.push(column);
    pivotRow += 1;
  }
  return { matrix, pivots };
}

function nullspaceFromRref(matrix, pivots) {
  const columns = matrix[0].length;
  const free = Array.from({ length: columns }, (_, j) => j)
    .filter((j) => !pivots.includes(j));
  return free.map((freeColumn) => {
    const vector = Array(columns).fill(0);
    vector[freeColumn] = 1;
    pivots.forEach((pivot, row) => {
      vector[pivot] = -matrix[row][freeColumn];
    });
    return vector;
  });
}

function fourSubspaces(A) {
  const reduced = rref(A);
  const reducedTranspose = rref(transpose(A));
  const rank = reduced.pivots.length;
  return {
    rowSpace: reduced.matrix.slice(0, rank),
    columnSpace: reduced.pivots.map(
      (column) => A.map((row) => row[column]),
    ),
    nullspace: nullspaceFromRref(reduced.matrix, reduced.pivots),
    leftNullspace: nullspaceFromRref(
      reducedTranspose.matrix,
      reducedTranspose.pivots,
    ),
    dimensions: {
      row: rank,
      column: rank,
      nullspace: A[0].length - rank,
      leftNullspace: A.length - rank,
    },
  };
}

const dot = (u, v) =>
  u.reduce((sum, value, i) => sum + value*v[i], 0);

function maxPairDot(first, second) {
  const values = first.flatMap((u) => second.map((v) => Math.abs(dot(u, v))));
  return values.length ? Math.max(...values) : 0;
}

const A = [
  [2, 4, 6, 4],
  [2, 5, 7, 6],
  [2, 3, 5, 2],
];
const spaces = fourSubspaces(A);
console.log({
  ...spaces,
  rowNullError: maxPairDot(spaces.rowSpace, spaces.nullspace),
  columnLeftError: maxPairDot(spaces.columnSpace, spaces.leftNullspace),
});

def transpose(A):
    return [list(column) for column in zip(*A)]

def rref(source, relative_tolerance=1e-12):
    matrix = [[float(value) for value in row] for row in source]
    row_scales = [max(abs(value) for value in row) for row in matrix]
    pivots = []
    pivot_row = 0
    for column in range(len(matrix[0])):
        if pivot_row == len(matrix):
            break
        best = max(
            range(pivot_row, len(matrix)),
            key=lambda row: (
                abs(matrix[row][column])/row_scales[row]
                if row_scales[row] else 0.0
            ),
        )
        best_ratio = (
            abs(matrix[best][column])/row_scales[best]
            if row_scales[best] else 0.0
        )
        if best_ratio <= relative_tolerance:
            continue
        matrix[pivot_row], matrix[best] = matrix[best], matrix[pivot_row]
        row_scales[pivot_row], row_scales[best] = (
            row_scales[best], row_scales[pivot_row]
        )
        scale = matrix[pivot_row][column]
        matrix[pivot_row] = [value/scale for value in matrix[pivot_row]]
        for row in range(len(matrix)):
            if row == pivot_row:
                continue
            factor = matrix[row][column]
            matrix[row] = [
                value - factor*pivot_value
                for value, pivot_value
                in zip(matrix[row], matrix[pivot_row])
            ]
        pivots.append(column)
        pivot_row += 1
    return matrix, pivots

def nullspace_from_rref(matrix, pivots):
    columns = len(matrix[0])
    free = [j for j in range(columns) if j not in pivots]
    basis = []
    for free_column in free:
        vector = [0.0]*columns
        vector[free_column] = 1.0
        for row, pivot in enumerate(pivots):
            vector[pivot] = -matrix[row][free_column]
        basis.append(vector)
    return basis

def four_subspaces(A):
    matrix, pivots = rref(A)
    matrix_t, pivots_t = rref(transpose(A))
    rank = len(pivots)
    return {
        "row_space": matrix[:rank],
        "column_space": [
            [row[column] for row in A] for column in pivots
        ],
        "nullspace": nullspace_from_rref(matrix, pivots),
        "left_nullspace": nullspace_from_rref(matrix_t, pivots_t),
        "dimensions": {
            "row": rank,
            "column": rank,
            "nullspace": len(A[0])-rank,
            "left_nullspace": len(A)-rank,
        },
    }

def dot(u, v):
    return sum(a*b for a, b in zip(u, v))

def max_pair_dot(first, second):
    values = [abs(dot(u, v)) for u in first for v in second]
    return max(values, default=0.0)

A = [
    [2, 4, 6, 4],
    [2, 5, 7, 6],
    [2, 3, 5, 2],
]
spaces = four_subspaces(A)
spaces["row_null_error"] = max_pair_dot(
    spaces["row_space"], spaces["nullspace"]
)
spaces["column_left_error"] = max_pair_dot(
    spaces["column_space"], spaces["left_nullspace"]
)
print(spaces)