本文主要是介绍c++ 线性代数 克·施密特（Gram Schmidt），希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！
克·施密特（Gram-Schmidt）正交化方法是一种将一组线性无关的向量转换为一组正交（垂直）向量的技术。该方法是线性代数中常用的工具，它的核心思想是将一组线性无关的向量集合通过减去它们在前面向量方向上的投影来得到一组正交的向量。
具体步骤如下：
- 给定一组线性无关的向量 {v₁, v₂, ..., vn}。
- 将第一个向量v₁单位化得到u₁： u₁ = v₁ / ||v₁||，其中||v₁||是v₁的模。
- 对于i从2到n： a. 计算第i个向量vᵢ与前i-1个单位化向量u₁, u₂, ..., uᵢ₋₁之间的投影： projᵢ = (vᵢ · u₁)u₁ + (vᵢ · u₂)u₂ + ... + (vᵢ · uᵢ₋₁)uᵢ₋₁ b. 计算第i个正交向量uᵢ： uᵢ = vᵢ - projᵢ c. 将uᵢ单位化： uᵢ = uᵢ / ||uᵢ||
通过这种方法，我们可以从一个线性无关的向量集合构造出一组正交的基向量，这在解决许多数学和工程问题时非常有用。Gram-Schmidt方法在解决线性方程组、矩阵分解、特征值问题等方面具有广泛的应用。
示例：
linearalgebra.h：
/* ----------------------- norm ----------------------- */
/* Given an array and its length, this function
computes the 2-norm of the array.
Input variables:
x : pointer to array for which the 2-norm should
be computed.
length: number of entries in x.
Features: This implementation has time complexity
O(length) and requires O(1) additional memory. */
double norm (double * x, int length) {
int i, length5;
double a, sum = 0;
length5 = length % 5;
for(i = 0; i < length5; i++) {
sum += x[i] * x[i];
}
for(; i < length; i += 5) {
sum += x[i] * x[i] + x[i + 1] * x[i + 1] + x[i + 2] * x[i + 2]
+ x[i + 3] * x[i + 3] + x[i + 4] * x[i + 4];
}
return sqrt(sum);
}
/* ----------------------- vec_copy ----------------------- */
/* Given two arrays of the same length and their length,
this function stores the values from the first array
in the second array.
Input variables:
x : pointer to array whose entries are to be
copied.
y : pointer to array in which the components
of x are to be stored.
length: number of entries in x and in y.
Features: This implementation has time complexity
O(length) and requires O(1) additional memory. */
void vec_copy (double * x, double * y, int length) {
int i, length5;
length5 = length % 5;
for(i = 0; i < length5; i++) {
y[i] = x[i];
}
for(; i < length; i += 5) {
y[i] = x[i];
y[i + 1] = x[i + 1];
y[i + 2] = x[i + 2];
y[i + 3] = x[i + 3];
y[i + 4] = x[i + 4];
}
}
/* ------------------- partialvec_copy ------------------- */
/* Given two arrays, the length of the second array, and
an index this function stores the values from the
subarray x[index : index + length] in the array
y[0 : length].
Input variables:
x : pointer to array whose entries are to be
copied.
y : pointer to array in which the components
of x are to be stored.
length: number of entries in y.
index : starting index of subarray of x to be
copied to y.
Example: Suppose x is a pointer to the array
{1, 2, 3, 4, 5}, y is a pointer to the array {0, 0, 0},
length = 3, and index = 2. Then after executing
partialvec_copy(x, y, 3, 2), the array pointed to by
y is now {3, 4, 5}.
Features: This implementation has time complexity
O(length) and requires O(1) additional memory. */
void partialvec_copy (double * x, double * y, int length, int index) {
int i, length5;
length5 = length % 5;
for(i = 0; i < length5; i++) {
y[i] = x[i + index];
}
for(; i < length; i += 5) {
y[i] = x[i + index];
y[i + 1] = x[i + index + 1];
y[i + 2] = x[i + index + 2];
y[i + 3] = x[i + index + 3];
y[i + 4] = x[i + index + 4];
}
}
/* ----------------------- scalar_div ----------------------- */
/* Given two arrays of the same length, their length, and a
scalar value this function divides the values from the
first array by the scalar value and stores the computed
number in the second array.
Input variables:
x : pointer to array whose components are to be
divided by r and stored in second array, y.
r : scalar used in division.
length: number of entries in x and in y.
y : pointer to array in which the components
of x are to be stored.
Features: This implementation has time complexity
O(length) and requires O(1) additional memory. */
void scalar_div (double * x, double r, int length, double * y) {
int i, length5;
length5 = length % 5;
for(i = 0; i < length5; i++) {
y[i] = x[i]/r;
}
for(; i < length; i += 5) {
y[i] = x[i]/r;
y[i + 1] = x[i + 1]/r;
y[i + 2] = x[i + 2]/r;
y[i + 3] = x[i + 3]/r;
y[i + 4] = x[i + 4]/r;
}
}
/* ----------------------- scalar_sub ----------------------- */
/* Given two arrays of the same length, their length, and a
scalar value this function multiplies the values from the
first array by the scalar value and then subtracts the
computed components from the components the second array.
Input variables:
x : pointer to array whose components are to be
multiplied by r then subtracted from the
components of the second array, y.
r : scalar used in multiplication.
length: number of entries in x and in y.
y : pointer to array in which the components
of x are to be stored.
Features: This implementation has time complexity
O(length) and requires O(1) additional memory. */
void scalar_sub (double * x, double r, int length, double * y) {
int i, length5;
length5 = length % 5;
for(i = 0; i < length5; i++) {
y[i] -= r * x[i];
}
for(; i < length; i += 5) {
y[i] -= r * x[i];
y[i + 1] -= r * x[i + 1];
y[i + 2] -= r * x[i + 2];
y[i + 3] -= r * x[i + 3];
y[i + 4] -= r * x[i + 4];
}
}
/* --------------------- partialscalar_sub --------------------- */
/* Given two arrays, the length of the second array, a scalar
value, and an index, this function multiplies the values
starting at the given index from the first array by the
scalar value and then subtracts the computed components from
the components the second array.
Input variables:
x : pointer to array whose components are to be
multiplied by r then subtracted from the
components of the second array, y.
r : scalar used in multiplication.
length: number of entries in y.
index :
y : pointer to array in which the components
of x are to be stored.
Example: Suppose x is a pointer to the array
{1, 2, 3, 4, 5}, y is a pointer to the array {0, 0, 0},
length = 3, r = -1, and index = 2. Then after executing
partialscalar_sub(x, -1, 3, 2, y), the array pointed to
by y is now {-3, -4, -5}.
Features: This implementation has time complexity
O(length) and requires O(1) additional memory. */
void partialscalar_sub (double * x, double r, int length,
int index, double * y)
{
int i, length5;
length5 = length % 5;
for(i = 0; i < length5; i++) {
y[i + index] -= r * x[i];
}
for(; i < length; i += 5) {
y[i + index] -= r * x[i];
y[i + index + 1] -= r * x[i + 1];
y[i + index + 2] -= r * x[i + 2];
y[i + index + 3] -= r * x[i + 3];
y[i + index + 4] -= r * x[i + 4];
}
}
/* --------------------- dot_product --------------------- */
/* Given two arrays of the same length and their length,
this function returns the dot product of the two
arrays.
Input variables:
x : pointer to first array.
y : pointer to second array.
length: number of entries in x and in y.
Features: This implementation has time complexity
O(length) and requires O(1) additional memory. */
double dot_product (double * x, double * y, int length) {
int i, length5;
double sum = 0;
length5 = length % 5;
for(i = 0; i < length5; i++) {
sum += x[i] * y[i];
}
for(; i < length; i += 5) {
sum += x[i] * y[i] + x[i + 1] * y[i + 1] + x[i + 2] * y[i + 2]
+ x[i + 3] * y[i + 3] + x[i + 4] * y[i + 4];
}
return sum;
}
/* ------------------ partialdot_product ------------------ */
/* Given two arrays of the same length, their length, and
an index this function returns the dot product of the
two subarrays x[index : length] and y[index : length].
Input variables:
x : pointer to first array.
y : pointer to second array.
length: number of entries in x and in y.
index : starting index for subarrays.
Example: Suppose x is a pointer to the array
{1, 2, 3, 4}, y is a pointer to the array {5, 6, 7, 8},
length = 4, and index = 2. Then the value returned by
executing partialdot_product(x, y, 4, 2) is 53, which
is computed by
x[2] * y[2] + x[3] * y[3] = 3 * 7 + 4 * 8
= 21 + 32
= 53.
Features: This implementation has time complexity
O(length) and requires O(1) additional memory. */
double partialdot_product (double * x, double * y, int length, int index) {
int i, length5;
double sum = 0;
length5 = length % 5;
for(i = index; i < length5; i++) {
sum += x[i] * y[i];
}
for(; i < length; i += 5) {
sum += x[i] * y[i] + x[i + 1] * y[i + 1] + x[i + 2] * y[i + 2]
+ x[i + 3] * y[i + 3] + x[i + 4] * y[i + 4];
}
return sum;
}
/* -------------------- subdot_product -------------------- */
/* Given two arrays, the length of the second array, and
an index this function returns the dot product of the
two subarrays x[index : index + length] and
y[0 : length]. It is necessary that index + length is
at most the length of the first array.
Input variables:
x : pointer to first array.
y : pointer to second array.
length: number of entries in y.
index : starting index for subarray of x.
Example: Suppose x is a pointer to the array
{1, 2, 3, 4, 5}, y is a pointer to the array
{-1, -2, -3}, length = 3, and index = 2. Then the value
returned by executing subdot_product(x, y, 3, 2) is 53,
which is computed by
x[2] * y[0] + x[3] * y[1] + x[4] * y[2]
= 3 * -1 + 4 * -2 + 5 * -3
= - 3 - 8 - 15
= -26.
Features: This implementation has time complexity
O(length) and requires O(1) additional memory. */
double subdot_product (double * x, double * y, int length, int index) {
int i, length5;
double sum = 0;
length5 = length % 5;
for(i = 0; i < length5; i++) {
sum += x[i + index] * y[i];
}
for(; i < length; i += 5) {
sum += x[i + index] * y[i] + x[i + index + 1] * y[i + 1]
+ x[i + index + 2] * y[i + 2]
+ x[i + index + 3] * y[i + 3]
+ x[i + index + 4] * y[i + 4];
}
return sum;
}
gramSchmidt.cpp：
#include <algorithm>
#include <cmath>
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include "linearalgebra.h"
/* ----------------------- gramSchmidt ----------------------- */
/* Given a matrix A of dimension m by n, this algorithm
computes a QR decomposition of A, where Q is a unitary
m by n matrix and R is a n by n upper triangular matrix
and A = QR.
Input variables:
a : pointer to array of arrays, the ith array of
which should correspond to the ith column of the
matrix A. During the algorithm, the columns of Q
will replace the columns of A.
r : pointer to array of arrays in which the ith
column of the upper triangular matrix R will be
stored in the ith subarray of r.
m : number of columns in A.
n : number of rows in A.
thin: TRUE => thin QR factorization computed
FALSE => full QR factorization computed
Features: This implementation has time complexity O(m n^2)
and requires O(1) additional memory.
Remarks: Due to the nature of the problem, if A is nearly
rank-deficient then the resulting columns of Q may not
exhibit the orthogonality property. */
void gramSchmidt (double ** a, double ** r, int m, int n, bool full) {
int i, j;
double anorm, tol = 10e-7;
for(i = 0; i < n; i++) {
r[i][i] = norm(a[i], m); // r_ii = ||a_i||
if(r[i][i] > tol) {
scalar_div(a[i], r[i][i], m, a[i]); // a_i = a_i/r_ii
}
else if(i == 0) { // set a[0] = [1 0 0 ... 0]^T
a[i][0] = 1;
for(j = 1; j < m; j++) {
a[i][j] = 0;
}
}
else{ // need to choose a_i orthogonal to < a_1, ... a_{i-1} >
for(j = 0; j < m; j++) {
a[i][j] = -a[0][i] * a[0][j];
}
a[i][i] += 1;
for(j = 1; j < i; j++) {
scalar_sub(a[j], a[j][i], m, a[i]);
}
anorm = norm(a[i], m);
scalar_div(a[i], anorm, m, a[i]);
}
for(j = i+1; j < n; j++) {
r[j][i] = dot_product(a[i], a[j], m); // r_ij = a_i*a_j
scalar_sub(a[i], r[j][i], m, a[j]); // a_j -= r_ij a_i
}
}
/* if full QR factorization requested, we choose remaining
columns of Q so that the m columns of Q form an
orthonormal set */
if(full) {
for(; i < m; i++) {
for(j = 0; j < m; j++) {
a[i][j] = -a[0][i] * a[0][j];
}
a[i][i] += 1;
for(j = 1; j < i; j++) {
scalar_sub(a[j], a[j][i], m, a[i]);
}
anorm = norm(a[i], m);
scalar_div(a[i], anorm, m, a[i]);
}
}
}
int main () {
int i, j, n, m, q_n, r_m;
bool full;
double x;
/* let user set the dimension of matrix A */
std::cout << "Enter the dimension m (where A is a m by n matrix): ";
std::cin >> m;
std::cout << "Enter the dimension n (where A is a m by n matrix): ";
std::cin >> n;
if(m != n) {
/* check if m < n */
if(m < n) {
printf("For a successful factorization, this implementation "
"requires n <= m.\nTerminating program.\n");
return 0;
}
/* let user choose either full or thin QR factorization */
std::cout << "Enter either 0 to compute a thin QR factorization"
<< std::endl;
std::cout << " or 1 to compute a full QR factorization: ";
std::cin >> full;
}
else { // else m == n so full and thin QR factorization are identical */
full = 1;
}
/* set dimensions of matrices Q and R based on full or thin QR */
if(full) { // Q is m by m and R is m by n
q_n = m;
r_m = m;
}
else { // Q is m by n and R is n by n
q_n = n;
r_m = n;
}
/* allocate memory for the matrices A and R */
double ** a = new double*[q_n];
double ** r = new double*[n];
for(i = 0; i < n; i++) {
a[i] = new double[m];
r[i] = new double[r_m];
}
for(; i < q_n; i++) {
a[i] = new double[m];
}
/* initialize the values in matrix A (only n columns regardless of
thin QR or full QR) */
for(i = 0; i < n; i++) {
for(j = i; j < m; j++) {
a[i][j] = j - i + 1; // this choice of values was arbitrary
}
}
/* print the matrix A before calling gramSchmidt */
std::cout << "A = " << std::endl;
for(i = 0; i < m; i++) {
for(j = 0; j < n; j++) {
printf("%9.6lg ", a[j][i]);
}
std::cout << std::endl;
}
std::cout << std::endl;
/* execute gramSchmidt to compute QR factorization */
gramSchmidt(a, r, m, n, full);
/* print the matrix Q resulting from gramSchmidt */
std::cout << "Q = " << std::endl;
for(i = 0; i < m; i++) {
for(j = 0; j < q_n; j++) {
if(a[j][i] >= 0) {
std::cout << " ";
}
printf("%9.6lg ", a[j][i]);
}
std::cout << std::endl;
}
std::cout << std::endl;
/* print the matrix R resulting from gramSchmidt */
std::cout << "R = " << std::endl;
for(i = 0; i < r_m; i++) {
for(j = 0; j < n; j++) {
printf("%9.6lg ", r[j][i]);
}
std::cout << std::endl;
}
std::cout << std::endl;
/* print numerical evidence that columns of Q are orthonormal */
printf("Numerical verification that {q_1, ..., q_%i} is an "
"orthonormal set:\n", q_n);
for(i = 0; i < q_n; i++) {
for(j = i; j < q_n; j++) {
x = dot_product(a[i], a[j], m);
printf("q_%i * q_%i = %lg\n", i + 1, j + 1, x);
}
}
/* free memory */
for(i = 0; i < n; i++) {
delete[] a[i];
delete[] r[i];
}
for(; i < q_n; i++) {
delete[] a[i];
}
delete[] a;
delete[] r;
return 0; // exit main
}
这篇关于c++ 线性代数 克·施密特（Gram Schmidt）的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
