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

e-maxx algorithm original: C/C++ #algorithm #emaxx #linear-algebra #math #matrix

Task text is translated from Russian for the selected interface language. Code is left unchanged.

Static

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

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

Algorithm

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

Implementation

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. it is required вычислить её определитель.

Algorithm

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

Implementation

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;
}

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

Vacancies for this task

Active vacancies with overlapping task tags are shown.

All vacancies

There are no active vacancies yet.