Kode C berikut adalah fungsi untuk melakukan inversi matrix singular dengan Metoda Marquardt untuk floating point matrix 2x2. Input: matrixA.txt dan matrixX.txt, dimana matrixA.txt adalah matrix singular atau ill-conditioned.
Di bagian bawah kode ini terdapat kode utama untuk perkalian matrix dan inversi marquardt dan membandingkan hasilnya.
Anda dapat melakukan test dimana:
matrixA.txt
2.0 2.0
3.0 3.0
matrixX.txt
0.5 0.5
0.5 0.5
#include<stdio.h>
void marquardt(float A[10][10],float X[10][10],float B[10][10],float lamda)
{
//compute AT=A'
float AT[10][10] ;
int i,j,k;
for(i=0;i<2;i++)
for(i=0;i<2;i++){
for(j=0;j<2;j++){
AT[j][i]=A[i][j];
}
}
//compute AAT=A'*A
float AAT[10][10];
float sum=0.0;
for(i=0;i<2;i++)
for(j=0;j<2;j++)
AAT[i][j]=0.0;
for(i=0;i<2;i++){
for(j=0;j<2;j++){
sum=0.0;
for(k=0;k<2;k++)
sum=sum+AT[i][k]*A[k][j];
AAT[i][j]=sum;
}
}
//create lamda*I
float LI[10][10];
for(i=0;i<2;i++)
{
for(j=0;j<2;j++)
{
if(i==j)
{
LI[i][j]=lamda;
}
else
{
LI[i][j]=0.0;
}
}
}
//compute AATLI=A'*A+lamda*I=AAT+LI
float AATLI[10][10];
for(i=0;i<2;i++)
{
for(j=0;j<2;j++)
{
AATLI[i][j] = (AAT[i][j] + LI[i][j]);
}
}
float ratio,a;
//compute inverse AATLI
{
for(i = 0; i < 2; i++){
for(j = 2; j < 2*2; j++){
if(i==(j-2))
AATLI[i][j] = 1.0;
else
AATLI[i][j] = 0.0;
}
}
for(i = 0; i < 2; i++){
for(j = 0; j < 2; j++){
if(i!=j){
ratio = AATLI[j][i]/AATLI[i][i];
for(k = 0; k < 2*2; k++){
AATLI[j][k] -= ratio * AATLI[i][k];
}
}
}
}
for(i = 0; i < 2; i++){
a = AATLI[i][i];
for(j = 2; j < 2*2; j++){
AATLI[i][j] /= a;
}
}
for(i = 0; i < 2; i++){
for(j = 0; j < 2; j++){
AATLI[i][j] = AATLI[i][2+j] ;
}
}
}
//compute C=inv(AATLI)*A'
float C[10][10];
for(i=0;i<2;i++)
for(j=0;j<2;j++)
C[i][j]=0.0;
for(i=0;i<2;i++){
for(j=0;j<2;j++){
sum=0.0;
for(k=0;k<2;k++)
sum=sum+AATLI[i][k]*AT[k][j];
C[i][j]=sum;
}
}
//compute B=inv(A'*A+lamda*I)*A'*X=C*X
//float B[10][10];
for(i=0;i<2;i++)
for(j=0;j<2;j++)
B[i][j]=0.0;
for(i=0;i<2;i++){
for(j=0;j<2;j++){
sum=0.0;
for(k=0;k<2;k++)
sum=sum+C[i][k]*X[k][j];
B[i][j]=sum;
}
}
}
//read matrix A
void readA(float A[10][10])
{
int i,j;
FILE *inputA=fopen("matrixA.txt","r");
for(i=0;i<2;++i){
for(j=0;j<2;++j){
fscanf(inputA,"%f",&A[i][j]);
}
}
fclose(inputA);
}
//read matrix X
void readX(float X[10][10])
{
int i,j;
FILE *inputX=fopen("matrixX.txt","r");
for(i=0;i<2;++i){
for(j=0;j<2;++j){
fscanf(inputX,"%f",&X[i][j]);
}
}
fclose(inputX);
}
//mult
void mult(float A[10][10],float B[10][10],float C[10][10])
{
int i,j,k;
float sum=0.0;
for(i=0;i<2;i++)
for(j=0;j<2;j++)
C[i][j]=0.0;
for(i=0;i<2;i++){
for(j=0;j<2;j++){
sum=0.0;
for(k=0;k<2;k++)
sum=sum+A[i][k]*B[k][j];
C[i][j]=sum;
}
}
}
//print
void print(float D[10][10])
{
int i,j;
for(i=0;i<2;++i)
{
for(j=0;j<2;++j)
{
printf("%f ",D[i][j]);
}
printf("\n");
}
}
main()
{
float A1[10][10],B1[10][10],C1[10][10],D1[10][10],X1[10][10];
int i,j,k;
readA(A1);//call matrix A
readX(X1);//call matrix X
mult(A1,X1,C1);//C=A*X
marquardt(A1,C1,D1,0.001);
print(X1); //print original input
print(D1); //print predicted input
}
No comments:
Post a Comment