Pada bagian ini saya telah menjelaskan teori SVD serta aplikasinya dengan Matlab serta SVDLIBC. Pada bagian ini saya akan menjelaskan penggunaan redsvd serta menghitung matrix inverse dengan input hasil SVD dengan menggunakan kode C.
Saya telah berhasil men-test program redsvd dan kode C berikut pada PC untuk matrix berukuran 5000x10000, sementara out of memory jika saya menggunakan Matlab.
Untuk dapat menggunakan redsvd anda harus memiliki platform Ubuntu 12.04.1 LTS. Setelah Ubuntu 12.04.1 LTS diinstall, anda perlu meng-install eigen3 melalui synaptic.
1. Install Synaptic
sudo apt-get update
sudo apt-get upgrade
sudo apt-get install synaptic
2. Install eigen3
Dash Home > search Synaptic > Synaptic Package Manager
Select libeigen3-dev dan libeigen3-doc lalu Apply
3. Install g++
sudo apt-get update
sudo apt-get install g++
4. Download redsvd di sini
Setelah didownload (pada kasus ini saya mendownload redsvd-0.2.0.tar.bz2)
tar xvjf redsvd-0.2.0.tar.bz2
cd redsvd-0.2.0
./waf configure
./waf
sudo ./waf install
5. Hitung SVD dengen redsvd
redsvd -i file1 -o file1 -r 10 -f dense
-i file1 artinya memanggil file1 sebagai input, -o file1 mengoutputkan file1.U, file-1.S, file1.V (V belum di-transpose) -r 10 artinya rank-10 atau mengoutputkan 10 non zero eigen terbesar -f dense artinya matrix input dense (bukan sparse). Dari pengalaman saya, nilai rank adalah nilai terbesar dari jumlah kolom atau baris dari matrix. Contoh matrix saya berukuran (5x4) saya pilih rank 5, matrix saya (3x7) saya pilih rank 7. Sebagai antisipasi, memilih rank yang besar tidak masalah sehingga tidak ada nilai eigen yang tidak di-outputkan.
Pada contoh ini saya memiliki file1 (matrix ini singular):
1.0 0.0 0.0
-2.0 0.0 0.0
4.0 6.0 1.0
6. Untuk memperoleh inverse dari matrix file1, copy-lah kode dibawah ini lalu beri nama inverse.c
7. Jalankan dengan perintah gcc -o inverse inverse.c
8. ./inverse
9. Perhatikan untuk mengubah nama file-file input, jumlah baris dan kolom dari input dan nama output.
10. Berikut adalah inverse dari matrix file1
~/Documents$ cat output.txt
0.200000 -0.400000 -0.000000
-0.129730 0.259460 0.162162
-0.021622 0.043243 0.027027
//This program calculates inverse of a matrix
//The input is U,S,V that are decomposed via SVD Decomposition
//Ainv=V*So*U' in which So=1/S
//by Agus Abdullah,PhD 2013.
#include<stdio.h>
#include <stdlib.h>
int main(void)
{
float **matrixU,*matrixS, **matrixV, **matrixI, **matrixD;
float **matrixSo,**matrixVSo, **matrixUT, **matrixVSoUT;
int row,col,i,j,nav1, nav2;
double eps;
//The matrix that you want to be inverted has row=??? and col=???
row=3;
col=3;
eps=0.0001;
//Read U
matrixU=(float **)malloc(row*sizeof(float));
for(i=0;i<row;i++)
matrixU[i]=(float *)malloc(row*sizeof(float));
FILE *inputU=fopen("file1.U","r");
for(i=0;i<row;i++)
{
for(j=0;j<row;j++)
{
fscanf(inputU,"%f",&matrixU[i][j]);
}
}
fclose(inputU);
printf("Step1 done... \n");
//Read V
matrixV=(float **)malloc(col*sizeof(float));
for(i=0;i<col;i++)
matrixV[i]=(float *)malloc(row*sizeof(float));
FILE *inputV=fopen("file1.V","r");
for(i=0;i<col;i++)
{
for(j=0;j<row;j++)
{
fscanf(inputV,"%f",&matrixV[i][j]);
}
}
fclose(inputV);
printf("Step2 done... \n");
//Read S
matrixS=(float *)malloc(row*sizeof(float));
FILE *inputS=fopen("file1.S","r");
for(i=0;i<row;i++)
{
fscanf(inputS,"%f",&matrixS[i]);
}
fclose(inputS);
//Matrix Identity
matrixI=(float **)malloc(row*sizeof(float));
for(i=0;i<row;i++)
matrixI[i]=(float *)malloc(row*sizeof(float));
for(i=0;i<row;i++)
{
for(j=0;j<row;j++)
{
if(i==j)
{
matrixI[i][j]=1.0;
}
else
{
matrixI[i][j]=0.0;
}
}
}
printf("Step3 done... \n");
//Matrix S Diagonalization
matrixD=(float **)malloc(row*sizeof(float));
for(i=0;i<row;i++)
matrixD[i]=(float *)malloc(row*sizeof(float));
for(i=0;i<row;i++)
{
for(j=0;j<row;j++)
{
matrixD[i][j]=matrixS[i]*matrixI[i][j];
}
}
printf("Step4 done... \n");
//Matrix 1/S
matrixSo=(float **)malloc(row*sizeof(float));
for(i=0;i<row;i++)
matrixSo[i]=(float *)malloc(row*sizeof(float));
for(i=0;i<row;i++)
{
for(j=0;j<row;j++)
if (matrixD[i][j] > eps)
{
matrixSo[i][j]=1/matrixD[i][j];
}
else
matrixSo[i][j]=0.0;
}
printf("Step5 done... \n");
//V*So --->V from redsvd is not transposed
matrixVSo=(float **)malloc(col*sizeof(float));
for(i=0;i<col;i++)
matrixVSo[i]=(float *)malloc(row*sizeof(float));
for(i=0;i<col;i++)
{
for(j=0;j<row;j++)
{ matrixVSo[i][j]=0;
for(nav1=0;nav1<row;nav1++)
matrixVSo[i][j]+=matrixV[i][nav1]*matrixSo[nav1][j];
}
}
printf("Step6 done... \n");
//Compute Transpose of U
matrixUT = malloc(row * sizeof(float));
for (i = 0; i < row; i++)
matrixUT[i] = malloc(row * sizeof(float));
for (i = 0; i < row; i++)
{
for (j = 0; j < row; j++)
matrixUT[j][i] = matrixU[i][j];
}
printf("Step7 done... \n");
//V*So*U'
matrixVSoUT=(float **)malloc(col*sizeof(float));
for(i=0;i<col;i++)
matrixVSoUT[i]=(float *)malloc(row*sizeof(float));
for(i=0;i<col;i++)
{
for(j=0;j<row;j++)
{ matrixVSoUT[i][j]=0;
for(nav2=0;nav2<row;nav2++)
matrixVSoUT[i][j]+=matrixVSo[i][nav2]*matrixUT[nav2][j];
}
}
printf("Step8 done... \n");
FILE *f;
f = fopen("output.txt","w");
for(i=0; i<col; i++)
{
for(j=0; j<row; j++)
{
fprintf(f,"\t%.4f",matrixVSoUT[i][j]);
}
fprintf(f,"\n");
}
fclose(f);
printf("completed! \n");
//Free Memory
for(i=0;i<row;i++)
free(matrixI[i]);
free(matrixI);
for(i=0;i<row;i++)
free(matrixU[i]);
free(matrixU);
for(i=0;i<col;i++)
free(matrixV[i]);
free(matrixV);
for(i=0;i<row;i++)
free(matrixD[i]);
free(matrixD);
free(matrixS);
for(i=0;i<row;i++)
free(matrixSo[i]);
free(matrixSo);
for(i=0;i<col;i++)
free(matrixVSo[i]);
free(matrixVSo);
for (i=0; i< row; i++)
free(matrixUT[i]);
free(matrixUT);
for(i=0;i<col;i++)
free(matrixVSoUT[i]);
free(matrixVSoUT);
}
No comments:
Post a Comment