#include "C
Matrix.h"
#include <algorithm>
#include <cmath>
#include <iomanip>
const double EPSILON = 1e-10; // 浮点数精度阈值
// 构造函数
C
Matrix::C
Matrix()
: rows(0), cols(0) {}
C
Matrix::C
Matrix(size_t r, size_t c, double initVal)
: rows(r), cols(c), data(r, std
::vector<double>(c, initVal)) {}
C
Matrix::C
Matrix(const std
::vector<std
::vector<double>>& d) {
if (d.empty() || d[0].empty()) throw std
::invalid_argument("Invalid
matrix dimensions");
rows = d.size();
cols = d[0].size();
data = d;
}
// 拷贝构造函数
C
Matrix::C
Matrix(const C
Matrix& other)
: rows(other.rows), cols(other.cols), data(other.data) {}
// 赋值运算符
C
Matrix& C
Matrix::operator=(const C
Matrix& other) {
if (this != &other) {
rows = other.rows;
cols = other.cols;
data = other.data;
}
return *this;
}
// 调整矩阵大小
void C
Matrix::resize(size_t newRows, size_t newCols, double initVal) {
data.resize(newRows);
for (auto& row
: data) {
row.resize(newCols, initVal);
}
rows = newRows;
cols = newCols;
}
// 矩阵加法
C
Matrix C
Matrix::operator+(const C
Matrix& rhs) const {
if (rows != rhs.rows || cols != rhs.cols) {
throw std
::invalid_argument("
Matrix dimensions do not match for addition");
}
C
Matrix result(rows, cols);
for (size_t i = 0; i < rows; ++i) {
for (size_t j = 0; j < cols; ++j) {
result.data[i][j] = data[i][j] + rhs.data[i][j];
}
}
return result;
}
// 矩阵减法
C
Matrix C
Matrix::operator-(const C
Matrix& rhs) const {
if (rows != rhs.rows || cols != rhs.cols) {
throw std
::invalid_argument("
Matrix dimensions do not match for subtraction");
}
C
Matrix result(rows, cols);
for (size_t i = 0; i < rows; ++i) {
for (size_t j = 0; j < cols; ++j) {
result.data[i][j] = data[i][j] - rhs.data[i][j];
}
}
return result;
}
// 矩阵乘法
C
Matrix C
Matrix::operator*(const C
Matrix& rhs) const {
if (cols != rhs.rows) {
throw std
::invalid_argument("
Matrix dimensions do not match for
multiplication");
}
C
Matrix result(rows, rhs.cols);
for (size_t i = 0; i < rows; ++i) {
for (size_t k = 0; k < cols; ++k) {
for (size_t j = 0; j < rhs.cols; ++j) {
result.data[i][j] += data[i][k] * rhs.data[k][j];
}
}
}
return result;
}
// 标量乘法（右乘）
C
Matrix C
Matrix::operator*(double scalar) const {
C
Matrix result(*this);
for (auto& row
: result.data) {
for (double& elem
: row) {
elem *= scalar;
}
}
return result;
}
// 标量除法
C
Matrix C
Matrix::operator/(double scalar) const {
if (std
::abs(scalar) < EPSILON) throw std
::invalid_argument("Division by zero");
return *this * (1.0 / scalar);
}
// 矩阵除法（乘以逆矩阵）
C
Matrix C
Matrix::operator/(const C
Matrix& rhs) const {
return *this * rhs.
inverse();
}
// 转置
C
Matrix C
Matrix::transpose() const {
C
Matrix result(cols, rows);
for (size_t i = 0; i < rows; ++i) {
for (size_t j = 0; j < cols; ++j) {
result.data[j][i] = data[i][j];
}
}
return result;
}
// 高斯消元法（全选主元）
void C
Matrix::gaussianElimination(std
::vector<std
::vector<double>>& mat, bool fullPivot) const {
size_t n = mat.size();
for (size_t i = 0; i < n; ++i) {
// 寻找主元
size_t maxRow = i, maxCol = i;
double maxVal = std
::abs(mat[i][i]);
if (fullPivot) {
for (size_t r = i; r < n; ++r) {
for (size_t c = i; c < n; ++c) {
if (std
::abs(mat[r][c]) > maxVal) {
maxVal = std
::abs(mat[r][c]);
maxRow = r;
maxCol = c;
}
}
}
// 交换行和列
if (maxRow != i) std
::swap(mat[i], mat[maxRow]);
if (maxCol != i) {
for (size_t r = 0; r < n; ++r) {
std
::swap(mat[r][i], mat[r][maxCol]);
}
}
}
// 消元
for (size_t j = i + 1; j < n; ++j) {
double factor = mat[j][i] / mat[i][i];
for (size_t k = i; k < mat[j].size(); ++k) {
mat[j][k] -= factor * mat[i][k];
}
}
}
}
// 求逆
C
Matrix C
Matrix::inverse() const {
if (rows != cols) throw std
::invalid_argument("Only square matrices can be inverted");
size_t n = rows;
std
::vector<std
::vector<double>> augMat(n, std
::vector<double>(2 * n, 0.0));
// 构造增广矩阵 [A | I]
for (size_t i = 0; i < n; ++i) {
for (size_t j = 0; j < n; ++j) {
augMat[i][j] = data[i][j];
}
augMat[i][n + i] = 1.0;
}
// 高斯-约当消元
gaussianElimination(augMat, true);
// 回代
for (int i = n - 1; i >= 0; --i) {
double diag = augMat[i][i];
if (std
::abs(diag) < EPSILON) throw std
::runtime_
error("
Matrix is
singular");
for (size_t j = 0; j < 2 * n; ++j) {
augMat[i][j] /= diag;
}
for (int k = i - 1; k >= 0; --k) {
double factor = augMat[k][i];
for (size_t j = 0; j < 2 * n; ++j) {
augMat[k][j] -= factor * augMat[i][j];
}
}
}
// 提取逆矩阵
C
Matrix inv(n, n);
for (size_t i = 0; i < n; ++i) {
for (size_t j = 0; j < n; ++j) {
inv.data[i][j] = augMat[i][j + n];
}
}
return inv;
}
// 求秩
size_t C
Matrix::rank() const {
std
::vector<std
::vector<double>> temp = data;
gaussianElimination(temp);
size_t rk = 0;
for (size_t i = 0; i < rows; ++i) {
bool nonZeroRow = false;
for (size_t j = 0; j < cols; ++j) {
if (std
::abs(temp[i][j]) > EPSILON) {
nonZeroRow = true;
break;
}
}
if (nonZeroRow) ++rk;
}
return rk;
}
// 标量乘法（左乘）
C
Matrix operator*(double scalar, const C
Matrix&
matrix) {
return
matrix * scalar;
}
// 输出运算符
std
::ostream& operator<<(std
::ostream& os, const C
Matrix&
matrix) {
for (size_t i = 0; i <
matrix.rows; ++i) {
for (size_t j = 0; j <
matrix.cols; ++j) {
os << std
::setw(10) <<
matrix.data[i][j] << " ";
}
os << std
::endl;
}
return os;
}
// 输入运算符
std
::istream& operator>>(std
::istream& is, C
Matrix&
matrix) {
size_t r, c;
is >> r >> c;
matrix.resize(r, c);
for (size_t i = 0; i < r; ++i) {
for (size_t j = 0; j < c; ++j) {
is >>
matrix.data[i][j];
}
}
return is;
}
简介缩短一点
本文探讨了在数据处理中遇到奇异矩阵问题的原因，包括解释变量过多或样本量不足等情况，以及如何通过调整样本量或减少变量数量来解决这一问题。
扫一扫
955
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
