我用 C 语言编写了一个程序，使用雅可比算法来计算矩阵的 SVD。我的程序中有 3 个函数:一个名为 jacobi_helper 的函数，它将执行 SVD；一个名为 jacobi 的函数，它在初始化几个数组后调用 jacobi_helper；以及一个调用 jacobi 的 main 函数。
我的问题在于，SVD 在 jacobi_helper 中似乎执行正确，但这些值似乎没有返回到 jacobi，即使它们是在 void 函数中传递的。最后，我尝试打印 SVD 的矩阵 U 和 V 以及奇异值 vector ，而我的困惑源于这样的事实:奇异值 vector 能够很好地打印，但矩阵却不能。
这是我的代码:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
double dot_product(double *a, double *b, int n){
double sum = 0;
for(int i = 0; i < n; i++){
sum = sum + a[i]*b[i];
}
return sum;
}
//Right matrix must be transpose of the matrix you want to multiply
void matrix_multiplication(left, right, output, n)
int n;
double left[n][n];
double right[n][n];
double output[n][n];
{
for(int i = 0; i < n; i++){
for(int j = 0; j < n; j++){
output[i][j] = dot_product(left[i], right[j], n);
}
}
}
void transpose(a, n, a_transpose)
int n;
double a[n][n];
double a_transpose[n][n];
{
for(int i = 0; i < n; i++){
for(int j = 0; j < n; j++){
a_transpose[i][j] = a[j][i];
}
}
}
void jacobi_helper(a, n, givens, givens_transpose, s, s_matrix, u, temp_u, u_output, v, temp_v, v_output)
int n; //Dimension of Matrix
double a[n][n]; //Input Matrix
double givens[n][n]; //Givens matrix to be formed in each step
double givens_transpose[n][n]; //Transpose of Givens Matrix
double s[n]; //Vector of singular values
double s_matrix[n][n]; //Matrix containing the singular values.
double u[n][n]; //U Matrix in SVD. (A = USV)
double temp_u[n][n]; //Temporary Matrix
double u_output[n][n];
double v[n][n];
double temp_v[n][n];
double v_output[n][n];
{
u = a;
//Set V to Identity Matrix
for(int i = 0; i < n; i++){
for(int j = 0; j < n; j++){
if(i == j){
v[i][j] = 1;
}
else{
v[i][j] = 0;
}
}
}
double n_squared = 0;
double value = 0;
double epsilon = 0.1;
for(int i = 0; i < n; i++){
for(int j = 0; j < n; j++){
n_squared = n_squared + (u[i][j]*u[i][j]);
}
}
do{
value = 0;
for(int i = 0; i < (n-1); i++){
for(int j = (i+1); j < n; j++){
double tmp1 = 0, tmp2 = 0, tmp4 = 0;
for(int k = 0; k < n; k++){
tmp1 = tmp1 + (u[k][i]*u[k][i]);
tmp2 = tmp2 + (u[k][i]*u[k][j]);
tmp4 = tmp4 + (u[k][j]*u[k][j]);
}
double tau = (tmp4 - tmp1)/(2*tmp2);
double t;
double t1;
double t2;
t1 = -tau + sqrt(1 + (tau*tau));
t2 = -tau - sqrt(1 + (tau*tau));
if(fabs(t1) < fabs(t2)){
t = t1;
}else{
t = t2;
}
double cos = 1/(sqrt(1+(t*t)));
double sin = cos*t;
//Make Givens Matrix
for(int k = 0; k < n; k++){
for(int l = 0; l < n; l++){
if(k == i && l == i){
givens[k][l] = cos;
}
else if (k == j && l == j){
givens[k][l] = cos;
}
else if ((k == i && l == j) && i < j){
givens[k][l] = -sin;
}
else if ((k == i && l == j) && i > j){
givens[k][l] = sin;
}
else if ((k == j && l == i) && i < j){
givens[k][l] = sin;
}
else if ((k == j && l == i) && i > j){
givens[k][l] = -sin;
}
else if(k == l){
givens[k][l] = 1;
}
else{
givens[k][l] = 0;
}
}
}
//Set U <- U*Givens
matrix_multiplication(u, givens, temp_u, n);
u = temp_u;
//Set V <- V*Givens
matrix_multiplication(v, givens, temp_v, n);
v = temp_v;
double temp = 0;
for(int k = 0; k < n; k++){
temp = temp + (u[k][i]*u[k][j]);
}
value = value + temp*temp;
}
}
}while(sqrt(value) >= (epsilon*epsilon)*(n_squared));
for(int i = 0; i < n; i++){
double sing_value = 0;
for(int k = 0; k < n; k++){
sing_value = sing_value + (u[k][i]*u[k][i]);
}
s[i] = sqrt(sing_value);
s_matrix[i][i] = 1/sqrt(sing_value);
}
matrix_multiplication(u, s_matrix, temp_u, n);
u = temp_u;
printf("%.16f %.16f %.16f %.16f\n",u[0][0],u[0][1],u[1][0],u[1][1]);
printf("%.16f %.16f %.16f %.16f\n",v[0][0],v[0][1],v[1][0],v[1][1]);
printf("%.16f %.16f\n\n", s[0], s[1]);
}
void jacobi(double *a, int n, double *s, double *u, double *v){
static double s_matrix[5000000];
static double givens[5000000];
static double givens_transpose[5000000];
static double temp_u[5000000];
static double temp_v[5000000];
static double u_output[5000000];
static double v_output[5000000];
jacobi_helper(a, n, givens, givens_transpose, s, s_matrix, u, temp_u, u_output, v, temp_v, v_output);
}
int main(int argc, const char * argv[]) {
static double a[4] = {2,2,-1,1};
int n = 2;
static double s[2];
static double u[4];
static double v[4];
jacobi(a, n, s, u, v);
printf("%.16f %.16f %.16f %.16f\n",u[0],u[1],u[2],u[3]);
printf("%.16f %.16f %.16f %.16f\n",v[0],v[1],v[2],v[3]);
printf("%.16f %.16f\n", s[0], s[1]);
return 0;
}
所以问题是，当我最后有 printf 语句时，s[0] 和 s[1] 在两种情况下都是相同的，但 u[0] 到 u[3] 和 v[0] 到v[3] 不一样。
从 jacobi_helper 函数内的 print 语句中我得到(正确的值):
u = 1.0000000000000000 0.0000000000000000 0.0000000000000000 1.0000000000000000
v = 0.7071067811865475 -0.7071067811865475 0.7071067811865475 0.7071067811865475
s = 2.8284271247461898 1.4142135623730949
但是从 main 函数我得到:
u = 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
v = 1.0000000000000000 0.0000000000000000 0.0000000000000000 1.0000000000000000
s = 2.8284271247461898 1.4142135623730949
谁能帮忙解决这个问题吗？如果这有用，我使用的算法是来自此网站的算法 6:http://www.cs.utexas.edu/users/inderjit/public_papers/HLA_SVD.pdf
最佳答案
jacobi_helper 使用参数 double *u 进行调用，该参数是一个指向可能存储结果的内存地址的指针，但 jacobi_helper 所做的第一件事是 u = a 的意思是:忘记内存地址 u 并使用 a 代替。您在指向参数的数组中设置了一些值v，但稍后您继续使用 v = temp_v。
数组变量只是指向内存地址的指针，为该数组变量赋值将用另一个地址替换该内存地址。如果您想更改数组变量指向的内存，则必须取消引用该地址。
代码
u[0] = 123;
u[1] = 345;
u[2] = 678;
更改u指向的内存中保存的数据。另一方面，
u = a;
只会更改 u，并且对 u[0] 的所有后续访问实际上将访问与 a[0] 相同的内存。
参数 double *u 的类型和参数 double u[n][n] 的类型也不匹配。
关于c - 数组在通过 void 函数后没有改变，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/28920070/