מרחב האפס

function rref(input, relativeTolerance=1e-12) {
  const matrix = input.map((row) => row.map(Number));
  const rows = matrix.length;
  const columns = matrix[0].length;
  const rowScales = matrix.map((row) => Math.max(...row.map(Math.abs)));
  const pivots = [];
  let pivotRow = 0;

  for (let column = 0; column < columns && pivotRow < rows; column += 1) {
    let best = -1;
    let bestRatio = 0;
    for (let row = pivotRow; row < rows; 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 < rows; 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 nullspaceBasis(A) {
  const { matrix, pivots } = rref(A);
  const columns = A[0].length;
  const free = Array.from(
    { length: columns },
    (_, index) => index,
  ).filter((index) => !pivots.includes(index));

  const basis = free.map((freeColumn) => {
    const vector = Array(columns).fill(0);
    vector[freeColumn] = 1;
    pivots.forEach((pivotColumn, row) => {
      vector[pivotColumn] = -matrix[row][freeColumn];
    });
    return vector;
  });
  return { pivots, free, basis };
}

function matVec(A, v) {
  return A.map((row) =>
    row.reduce((sum, value, j) => sum + value*v[j], 0),
  );
}

const A = [
  [2, 4, 6, 4],
  [2, 5, 7, 6],
  [2, 3, 5, 2],
];
const result = nullspaceBasis(A);
console.log({
  pivotVariables: result.pivots.map((i) => i + 1),
  freeVariables: result.free.map((i) => i + 1),
  basis: result.basis,
  verification: result.basis.map((vector) => matVec(A, vector)),
});

def rref(source, relative_tolerance=1e-12):
    matrix = [[float(value) for value in row] for row in source]
    rows, columns = len(matrix), len(matrix[0])
    row_scales = [max(abs(value) for value in row) for row in matrix]
    pivots = []
    pivot_row = 0

    for column in range(columns):
        if pivot_row == rows:
            break
        best = max(
            range(pivot_row, rows),
            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(rows):
            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_basis(A):
    matrix, pivots = rref(A)
    columns = len(A[0])
    free = [index for index in range(columns) if index not in pivots]
    basis = []
    for free_column in free:
        vector = [0.0]*columns
        vector[free_column] = 1.0
        for row, pivot_column in enumerate(pivots):
            vector[pivot_column] = -matrix[row][free_column]
        basis.append(vector)
    return {"pivots": pivots, "free": free, "basis": basis}

def mat_vec(A, vector):
    return [
        sum(value*vector[j] for j, value in enumerate(row))
        for row in A
    ]

A = [
    [2, 4, 6, 4],
    [2, 5, 7, 6],
    [2, 3, 5, 2],
]
result = nullspace_basis(A)
print({
    "pivot_variables": [i+1 for i in result["pivots"]],
    "free_variables": [i+1 for i in result["free"]],
    "basis": result["basis"],
    "verification": [mat_vec(A, vector) for vector in result["basis"]],
})