← Static tasks

E126. Нахождение ранга матрицы

emaxx algorithm

#algorithm#emaxx#linear-algebra#math#matrix

Task

Источник: e-maxx.ru/algo, страница PDF 415.

Ранг матрицы - это наибольшее число линейно независимых строк/столбцов матрицы. Ранг определён не только для квадратных матриц; пусть матрица прямоугольна и имеет размер NxM. Также ранг матрицы можно определить как наибольший из порядков миноров матрицы, отличных от нуля. Заметим, что если матрица квадратная и её определитель отличен от нуля, то ранг равен N(=M), иначе он будет меньше. В общем случае, ранг матрицы не превосходит min(N,M).

Алгоритм

Искать ранг можно с помощью модифицированного метода Гаусса. Будем выполнять абсолютно те же самые операции, что и при решении системы или нахождении её определителя, но если на каком-либо шаге в i-ом столбце среди невыбранных до этого строк нет ненулевых, то мы этот шаг пропускаем, а ранг уменьшаем на единицу (изначально ранг полагаем равным max(N,M)). Иначе, если мы нашли на i-ом шаге строку с ненулевым элементом в i- ом столбце, то помечаем эту строку как выбранную, и выполняем обычные операции отнимания этой строки от остальных.

Реализация

const double EPS = 1E-9;

int rank = max(n,m);

vector<char> line_used (n);

for (int i=0; i<m; ++i) {
int j;
for (j=0; j<n; ++j)
if (!line_used[j] && abs(a[j][i]) > EPS)

break;

if (j == n)

--rank;

else {

line_used[j] = true;

for (int p=i+1; p<m; ++p)

a[j][p] /= a[j][i];

for (int k=0; k<n; ++k)
if (k != j && abs (a[k][i]) > EPS)
for (int p=i+1; p<m; ++p)

a[k][p] -= a[j][p] * a[k][i];

} }

Вычисление определителя матрицы методом Гаусса

Пусть дана квадратная матрица A размером NxN. Требуется вычислить её определитель.

Алгоритм

Воспользуемся идеями метода Гаусса решения систем линейных уравнений. Будем выполнять те же самые действия, что и при решении системы линейных уравнений, исключив только деление текущей строки на a[i][i] (точнее, само деление можно выполнять, но подразумевая, что число выносится за знак определителя). Тогда все операции, которые мы будем производить с матрицей, не будут изменять величину определителя матрицы, за исключением, быть может, знака (мы только обмениваем местами две строки, что меняет знак на противоположный, или прибавляем одну строку к другой, что не меняет величину определителя). Но матрица, к которой мы приходим после выполнения алгоритма Гаусса, является диагональной, и определитель её равен произведению элементов, стоящих на диагонали. Знак, как уже говорилось, будет определяться количеством обменов строк (если их нечётное, то знак определителя следует изменить на противоположный). Таким образом, мы можем с помощью алгоритма Гаусса вычислять определитель матрицы за O (N3). Осталось только заметить, что если в какой-то момент мы не найдём в текущем столбце ненулевого элемента, то алгоритм следует остановить и вернуть 0.

Реализация

const double EPS = 1E-9;

int n;

vector < vector<double> > a (n, vector<double> (n));

... чтение n и a ...

double det = 1;

for (int i=0; i<n; ++i) {
int k = i;
for (int j=i+1; j<n; ++j)
if (abs (a[j][i]) > abs (a[k][i]))

k = j;

if (abs (a[k][i]) < EPS) {

det = 0;

break;

}

swap (a[i], a[k]);

if (i != k)

det = -det;

det *= a[i][i];

for (int j=i+1; j<n; ++j)

a[i][j] /= a[i][i];

for (int j=0; j<n; ++j)
if (j != i && abs (a[j][i]) > EPS)
for (int k=i+1; k<n; ++k)

a[j][k] -= a[i][k] * a[j][i];

}

cout << det;

C# solution

auto-draft, review before submit
using System;
using System.Collections.Generic;
using System.Linq;

public static class AlgorithmDraft
{
    // Auto-generated C# draft from the original e-maxx C/C++ listing. Review before production use.
    const double EPS = 1E-9;
    int rank = max(n,m);
    List<char> line_used (n);
    for (int i=0; i<m; ++i) {
            int j;
            for (j=0; j<n; ++j)
                    if (!line_used[j] && abs(a[j][i]) > EPS)
                            break;
            if (j == n)
                    --rank;
            else {
                    line_used[j] = true;
                    for (int p=i+1; p<m; ++p)
                            a[j][p] /= a[j][i];
                    for (int k=0; k<n; ++k)
                            if (k != j && abs (a[k][i]) > EPS)
                                    for (int p=i+1; p<m; ++p)
                                            a[k][p] -= a[j][p] * a[k][i];
            }
    }
    const double EPS = 1E-9;
    int n;
    vector < vector<double> > a (n, vector<double> (n));
    ... чтение n и a ...
    double det = 1;
    for (int i=0; i<n; ++i) {
            int k = i;
            for (int j=i+1; j<n; ++j)
                    if (abs (a[j][i]) > abs (a[k][i]))
                            k = j;
            if (abs (a[k][i]) < EPS) {
                    det = 0;
                    break;
            }
            swap (a[i], a[k]);
            if (i != k)
                    det = -det;
            det *= a[i][i];
            for (int j=i+1; j<n; ++j)
                    a[i][j] /= a[i][i];
            for (int j=0; j<n; ++j)
                    if (j != i && abs (a[j][i]) > EPS)
                            for (int k=i+1; k<n; ++k)
                                    a[j][k] -= a[i][k] * a[j][i];
    }
    Console.WriteLine( det;
}

C++ solution

matched/original
const double EPS = 1E-9;
int rank = max(n,m);
vector<char> line_used (n);
for (int i=0; i<m; ++i) {
        int j;
        for (j=0; j<n; ++j)
                if (!line_used[j] && abs(a[j][i]) > EPS)
                        break;
        if (j == n)
                --rank;
        else {
                line_used[j] = true;
                for (int p=i+1; p<m; ++p)
                        a[j][p] /= a[j][i];
                for (int k=0; k<n; ++k)
                        if (k != j && abs (a[k][i]) > EPS)
                                for (int p=i+1; p<m; ++p)
                                        a[k][p] -= a[j][p] * a[k][i];
        }
}
const double EPS = 1E-9;
int n;
vector < vector<double> > a (n, vector<double> (n));
... чтение n и a ...
double det = 1;
for (int i=0; i<n; ++i) {
        int k = i;
        for (int j=i+1; j<n; ++j)
                if (abs (a[j][i]) > abs (a[k][i]))
                        k = j;
        if (abs (a[k][i]) < EPS) {
                det = 0;
                break;
        }
        swap (a[i], a[k]);
        if (i != k)
                det = -det;
        det *= a[i][i];
        for (int j=i+1; j<n; ++j)
                a[i][j] /= a[i][i];
        for (int j=0; j<n; ++j)
                if (j != i && abs (a[j][i]) > EPS)
                        for (int k=i+1; k<n; ++k)
                                a[j][k] -= a[i][k] * a[j][i];
}
cout << det;

Java solution

auto-draft, review before submit
import java.util.*;
import java.math.*;

public class AlgorithmDraft {
    // Auto-generated Java draft from the original e-maxx C/C++ listing. Review before production use.
    const double EPS = 1E-9;
    int rank = max(n,m);
    ArrayList<Character> line_used (n);
    for (int i=0; i<m; ++i) {
            int j;
            for (j=0; j<n; ++j)
                    if (!line_used[j] && abs(a[j][i]) > EPS)
                            break;
            if (j == n)
                    --rank;
            else {
                    line_used[j] = true;
                    for (int p=i+1; p<m; ++p)
                            a[j][p] /= a[j][i];
                    for (int k=0; k<n; ++k)
                            if (k != j && abs (a[k][i]) > EPS)
                                    for (int p=i+1; p<m; ++p)
                                            a[k][p] -= a[j][p] * a[k][i];
            }
    }
    const double EPS = 1E-9;
    int n;
    vector < vector<double> > a (n, vector<double> (n));
    ... чтение n и a ...
    double det = 1;
    for (int i=0; i<n; ++i) {
            int k = i;
            for (int j=i+1; j<n; ++j)
                    if (abs (a[j][i]) > abs (a[k][i]))
                            k = j;
            if (abs (a[k][i]) < EPS) {
                    det = 0;
                    break;
            }
            swap (a[i], a[k]);
            if (i != k)
                    det = -det;
            det *= a[i][i];
            for (int j=i+1; j<n; ++j)
                    a[i][j] /= a[i][i];
            for (int j=0; j<n; ++j)
                    if (j != i && abs (a[j][i]) > EPS)
                            for (int k=i+1; k<n; ++k)
                                    a[j][k] -= a[i][k] * a[j][i];
    }
    System.out.println( det;
}

Explanation

Материал разбит как алгоритмическая задача: изучить постановку, понять асимптотику и реализовать алгоритм на выбранном языке.