RcdMathLib_doc
Open Source Library for Linear and Non-linear Algebra
moore_penrose_pseudo_inverse.c
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
3  * 2020 Freie Universität Berlin
4  *
5  * This file is subject to the terms and conditions of the GNU Lesser General
6  * Public License v2.1. See the file LICENSE in the top level directory for more
7  * details.
8  */
9 
24 #include <stdbool.h>
25 #include <stdio.h>
26 
28 #include "matrix.h"
29 #include "svd.h"
30 
31 static int32_t min(int32_t a, int32_t b);
32 static int8_t moore_penrose_pinv(uint8_t m, uint8_t n, matrix_t A[m][n],
33  uint8_t u_m, uint8_t u_n, matrix_t U[u_m][u_n],
34  matrix_t S[u_n][n], matrix_t V[n][n],
35  uint8_t s_length, matrix_t s[],
36  matrix_t pinv_A[n][m]);
37 
38 int8_t moore_penrose_get_pinv(uint8_t m, uint8_t n, matrix_t A[m][n],
39  matrix_t pinv_A[n][m])
40 {
41  matrix_dim_t udim;
42  int8_t answer = 0;
43 
44  if (m < n) { /*
45  * first compute the pseudo inverse of the transposed of the
46  * matrix A, than transpose the computed pseudo inverse.
47  */
48  uint8_t trans_m = n;
49  uint8_t trans_n = m;
50  matrix_t trans_A[n][m];
51  matrix_t trans_A_pinv[trans_n][trans_m];
52 
53  matrix_transpose(m, n, A, trans_A);
54  svd_get_U_dim(trans_m, trans_n, &udim);
55  matrix_t U[udim.row_num][udim.col_num];
56 
57  matrix_t S[udim.col_num][trans_n];
58  matrix_t V[trans_n][trans_n];
59  uint8_t s_length = svd_get_single_values_num(trans_m, trans_n);
60  matrix_t s[s_length];
61 
62  answer = moore_penrose_pinv(trans_m, trans_n, trans_A,
63  udim.row_num, udim.col_num, U, S, V,
64  s_length, s, trans_A_pinv);
65  matrix_transpose(trans_n, trans_m, trans_A_pinv, pinv_A);
66 
67  return answer;
68  }
69  else { /* compute the pseudo inverse */
70  svd_get_U_dim(m, n, &udim);
71  matrix_t U[udim.row_num][udim.col_num];
72  matrix_t S[udim.col_num][n];
73  matrix_t V[n][n];
74  uint8_t s_length = svd_get_single_values_num(m, n);
75  matrix_t s[s_length];
76  answer = moore_penrose_pinv(m, n, A,
77  udim.row_num, udim.col_num, U, S, V,
78  s_length, s, pinv_A);
79 
80  return answer;
81  }
82 
83 }
84 
101 static int8_t moore_penrose_pinv(uint8_t m, uint8_t n, matrix_t A[m][n],
102  uint8_t u_m, uint8_t u_n, matrix_t U[u_m][u_n],
103  matrix_t S[u_n][n], matrix_t V[n][n],
104  uint8_t s_length, matrix_t s[],
105  matrix_t pinv_A[n][m])
106 {
107 
108  matrix_dim_t u_dim;
109  matrix_t rec_s[s_length];
110 
111  uint8_t i, j, k, col_min;
112 
113  if ((m > MAX_ROW_NUM) | (n > MAX_COL_NUM)) {
114 
115  printf(
116  "The maximal, allowed number of the rows and columns is %d\n",
117  MAX_ROW_NUM);
119  }
120 
121  if (m < n) {
122 
123  puts("Please, give the transposed matrix");
125  }
126  else {
127  svd_get_U_dim(m, n, &u_dim);
128  s_length = svd_get_single_values_num(m, n);
129 
130  svd(m, n, A, u_dim.row_num, u_dim.col_num, U, S, V, s_length,
131  s);
132  }
133 
134  if (matrix_get_rank(m, n, s, s_length) < 1) {
135  puts("The minimum rank-value of a matrix is one");
137  }
138 
139  svd_get_reciproc_singular_values(m, n, s_length, s, rec_s);
140 
141  col_min = min(n, u_dim.col_num);
142 
143  matrix_clear(n, m, pinv_A);
144 
145  for (i = 0; i < n; i++) {
146  for (j = 0; j < u_dim.row_num; j++) {
147  for (k = 0; k < col_min; k++) {
148  pinv_A[i][j] += V[i][k] * rec_s[k] * U[j][k];
149 
150  }
151  }
152  }
153 
155 }
156 
157 static int32_t min(int32_t a, int32_t b)
158 {
159  if (a < b) {
160  return a;
161  }
162  else {
163  return b;
164  }
165 }
166 
167 void moore_penrose_pinv_compute_print(uint8_t m, uint8_t n,
168  matrix_t matrix[m][n], uint8_t i)
169 {
170  matrix_t pinv_mat[n][m];
171 
172  moore_penrose_get_pinv(m, n, matrix, pinv_mat);
173  printf("pinv%d = ", i);
174  matrix_flex_print(n, m, pinv_mat, 11, 7);
175  puts("");
176 }
MAX_ROW_NUM
#define MAX_ROW_NUM
The maximal row number allowed.
Definition: moore_penrose_pseudo_inverse.h:31
MOORE_PENROSE_INVALID_RANK_VALUE
#define MOORE_PENROSE_INVALID_RANK_VALUE
Invalid rank value of a matrix.
Definition: moore_penrose_pseudo_inverse.h:57
matrix_clear
void matrix_clear(uint8_t m, uint8_t n, matrix_t matrix[m][n])
Clear all the elements of the vector.
Definition: matrix.c:46
matrix_dim_t::row_num
uint8_t row_num
the row number
Definition: matrix.h:61
svd
void svd(uint8_t m, uint8_t n, matrix_t A[m][n], uint8_t u_m, uint8_t u_n, matrix_t U[u_m][u_n], matrix_t S[u_n][n], matrix_t V[n][n], uint8_t sing_vec_length, matrix_t singl_values_vec[])
Compute the Singular-Value Decomposition (SVD) of a matrix.
Definition: svd.c:119
MAX_COL_NUM
#define MAX_COL_NUM
The maximal column number allowed.
Definition: moore_penrose_pseudo_inverse.h:36
matrix_transpose
void matrix_transpose(uint8_t m, uint8_t n, matrix_t src_matrix[m][n], matrix_t dest_matrix[n][m])
Computes the transpose of a matrix.
Definition: matrix.c:56
MOORE_PENROSE_PSEUDO_GIVE_MATRIX_TRANSPOSE
#define MOORE_PENROSE_PSEUDO_GIVE_MATRIX_TRANSPOSE
The transposed matrix should be delivered.
Definition: moore_penrose_pseudo_inverse.h:52
MOORE_PENROSE_PSEUDO_COMP_SUCCESS
#define MOORE_PENROSE_PSEUDO_COMP_SUCCESS
The Moore–Penrose inverse is successfully completed.
Definition: moore_penrose_pseudo_inverse.h:42
matrix_get_rank
uint8_t matrix_get_rank(uint8_t m, uint8_t n, matrix_t singl_values_arr[], uint8_t length)
Compute the rank of a matrix.
Definition: matrix.c:312
matrix.h
Matrix computations.
matrix_flex_print
void matrix_flex_print(uint8_t m, uint8_t n, matrix_t matrix[m][n], uint8_t before_dec, uint8_t after_dec)
Display the values of the matrix elements.
Definition: matrix.c:220
moore_penrose_pseudo_inverse.h
Moore–Penrose algorithm to compute the pseudo-inverse of a matrix.
matrix_dim_t
A structure to define the row and column number of a matrix.
Definition: matrix.h:60
svd.h
Algorithm for the Singular Value Decomposition (SVD).
svd_get_single_values_num
uint8_t svd_get_single_values_num(uint8_t m, uint8_t n)
Calculate the number of the singular values.
Definition: svd.c:114
svd_get_reciproc_singular_values
void svd_get_reciproc_singular_values(uint8_t m, uint8_t n, uint8_t length, matrix_t singl_values_arr[], matrix_t recip_singl_values_arr[])
Compute the reciprocal singular values.
Definition: svd.c:747
matrix_t
#define matrix_t
Define the data type of the matrix elements.
Definition: matrix.h:38
MOORE_PENROSE_PSEUDO_MAX_ALLOW_ROW_COL_EXCEEED
#define MOORE_PENROSE_PSEUDO_MAX_ALLOW_ROW_COL_EXCEEED
The maximal row number allowed is exceeded.
Definition: moore_penrose_pseudo_inverse.h:47
moore_penrose_pinv_compute_print
void moore_penrose_pinv_compute_print(uint8_t m, uint8_t n, matrix_t matrix[m][n], uint8_t i)
Compute and print the Moore–Penrose inverse of a matrix.
Definition: moore_penrose_pseudo_inverse.c:167
moore_penrose_get_pinv
int8_t moore_penrose_get_pinv(uint8_t m, uint8_t n, matrix_t A[m][n], matrix_t pinv_A[n][m])
Calculate the Moore–Penrose inverse of a rectangular matrix.
Definition: moore_penrose_pseudo_inverse.c:38
svd_get_U_dim
void svd_get_U_dim(uint8_t m, uint8_t n, matrix_dim_t *u_dim)
Calculate the dimension of the matrix U.
Definition: svd.c:96
matrix_dim_t::col_num
uint8_t col_num
the column number
Definition: matrix.h:62