ベクトルと行列

永井 忠一 2020.5.30

ベクトル

ベクトルの演算

ベクトル同士の和と差。スカラーとの積（スカラー倍）

C 言語 Fortran Maxima Octave Python NumPy R 言語
OpenGL Mathematics (GLM) BLAS 
#include <stdio.h>
#include <glm/vec3.hpp>

int main(int argc, char *argv[]) {
	using namespace glm;
	vec3 u = vec3(1.f, 2.f, 3.f);
	vec3 v = vec3(4.f, 5.f, 6.f);

	vec3 ans = u + v;
	printf("%f %f %f\n", ans[0], ans[1], ans[2]);

	ans = u - v;
	printf("%f %f %f\n", ans[0], ans[1], ans[2]);

	ans = vec3(2.f)*u;
	printf("%f %f %f\n", ans[0], ans[1], ans[2]);

	return 0;
}
 
#include <stdio.h>
#include <cblas.h>

int main(int argc, char *argv[]) {
#define N 3
	{
		double alpha = 1.;
		double x[N] = { 1., 2., 3. };
		double y[N] = { 4., 5., 6. };
		cblas_daxpy(N, alpha, x, 1, y, 1);
		printf("%f %f %f\n", y[0], y[1], y[2]);
	}
	{
		double alpha = 1.;
		double x[N] = { 1., 2., 3. };
		double y[N] = { 4., 5., 6. };
		y[0] = -y[0]; y[1] = -y[1]; y[2] = -y[2];
		cblas_daxpy(N, alpha, x, 1, y, 1);
		printf("%f %f %f\n", y[0], y[1], y[2]);
	}
	{
		double alpha = 2.;
		double x[N] = { 1., 2., 3. };
		double y[N] = { 0., 0., 0. };
		cblas_daxpy(N, alpha, x, 1, y, 1);
#undef N
		printf("%f %f %f\n", y[0], y[1], y[2]);
	}
	return 0;
}
 
PROGRAM vec_prog1
  IMPLICIT NONE
  REAL(8), DIMENSION(3) :: u = (/ 1.D0, 2.D0, 3.D0 /)
  REAL(8), DIMENSION(3) :: v = (/ 4.D0, 5.D0, 6.D0 /)
  REAL(8), DIMENSION(3) :: ans

  ans = u + v
  WRITE(*,*) ans

  ans = u - v
  WRITE(*,*) ans

  ans = 2.D0*u
  WRITE(*,*) ans
END PROGRAM vec_prog1
 
(%i1) u: [1, 2, 3];
(%o1)                              [1, 2, 3]
(%i2) v: [4, 5, 6];
(%o2)                              [4, 5, 6]
(%i3) u + v;
(%o3)                              [5, 7, 9]
(%i4) u - v;
(%o4)                           [- 3, - 3, - 3]
(%i5) 2*u;
(%o5)                              [2, 4, 6]

octave:1> u = [1 2 3]
u =

   1   2   3

octave:2> v = [4 5 6]
v =

   4   5   6

octave:3> u + v
ans =

   5   7   9

octave:4> u - v
ans =

  -3  -3  -3

octave:5> 2*u
ans =

   2   4   6


In [1]: import numpy as np

In [2]: u = np.array([1, 2, 3])

In [3]: v = np.array([4, 5, 6])

In [4]: u + v
Out[4]: array([5, 7, 9])

In [5]: u - v
Out[5]: array([-3, -3, -3])

In [6]: 2*u
Out[6]: array([2, 4, 6])


 > u <- c(1, 2, 3)
 > v <- c(4, 5, 6)
 > u + v
 [1] 5 7 9
 > u - v
 [1] -3 -3 -3
 > 2*u
 [1] 2 4 6

$ g++ vec_prog1.cpp && ./a.out
5.000000 7.000000 9.000000
-3.000000 -3.000000 -3.000000
2.000000 4.000000 6.000000

$ clang vec_prog1.cpp && ./a.out
5.000000 7.000000 9.000000
-3.000000 -3.000000 -3.000000
2.000000 4.000000 6.000000
 
$ g++ la_prog1.cpp -lblas && ./a.out
5.000000 7.000000 9.000000
-3.000000 -3.000000 -3.000000
2.000000 4.000000 6.000000
 
$ gfortran vec_prog1.f90 && ./a.out
   5.0000000000000000        7.0000000000000000        9.0000000000000000     
  -3.0000000000000000       -3.0000000000000000       -3.0000000000000000     
   2.0000000000000000        4.0000000000000000        6.0000000000000000

ベクトルの内積(dot product)

C 言語 Fortran Maxima Octave Python R 言語
GLM BLAS 
#include <stdio.h>
#include <glm/glm.hpp>

int main(int argc, char *argv[]) {
	using namespace glm;
	vec3 u = vec3(1.f, 2.f, 3.f);
	vec3 v = vec3(4.f, 5.f, 6.f);

	float ans = dot(u, v);
	printf("%f\n", ans);

	return 0;
}
 
#include <stdio.h>
#include <cblas.h>

int main(int argc, char *argv[]) {
#define N 3
	double x[N] = { 1., 2., 3. };
	double y[N] = { 4., 5., 6. };
	double dot = cblas_ddot(N, x, 1, y, 1);
#undef N
	printf("%f\n", dot);

	return 0;
}
 
PROGRAM vec_prog2
  IMPLICIT NONE
  REAL(8), DIMENSION(3) :: u = (/ 1.D0, 2.D0, 3.D0 /)
  REAL(8), DIMENSION(3) :: v = (/ 4.D0, 5.D0, 6.D0 /)
  REAL(8) :: ans

  ans = DOT_PRODUCT(u, v)
  WRITE(*,*) ans
END PROGRAM vec_prog2
 
(%i1) u: [1, 2, 3];
(%o1)                              [1, 2, 3]
(%i2) v: [4, 5, 6];
(%o2)                              [4, 5, 6]
(%i3) u.v;
(%o3)                                 32
 
octave:1> u = [1 2 3]
u =

   1   2   3

octave:2> v = [4 5 6]
v =

   4   5   6

octave:3> dot(u, v)
ans =  32
octave:4> u*v'
ans =  32
 
In [1]: import numpy as np

In [2]: u = np.array([1, 2, 3])

In [3]: v = np.array([4, 5, 6])

In [4]: u.dot(v)
Out[4]: 32

In [5]: np.inner(u, v)
Out[5]: 32


 > u <- c(1, 2, 3)
 > v <- c(4, 5, 6)
 > u%*%v
      [,1]
 [1,]   32
 
$ g++ vec_prog2.cpp && ./a.out
32.000000

$ clang vec_prog2.cpp && ./a.out
32.000000
 
$ g++ la_prog2.cpp -lblas && ./a.out
32.000000
 
$ gfortran vec_prog2.f90 && ./a.out
   32.000000000000000

ベクトルの外積(cross product)

C言語 GLM Maxima Octave Python R言語 
#include <stdio.h>
#include <glm/glm.hpp>

int main(int argc, char *argv[]) {
	using namespace glm;
	vec3 u = vec3(1.f, 2.f, 3.f);
	vec3 v = vec3(4.f, 5.f, 6.f);

	vec3 ans = cross(u, v);
	printf("%f %f %f\n", ans[0], ans[1], ans[2]);

	return 0;
}
 
(%i1) load(vect);
(%o1)           /usr/share/maxima/5.37.2/share/vector/vect.mac
(%i2) u: [1, 2, 3];
(%o2)                              [1, 2, 3]
(%i3) v: [4, 5, 6];
(%o3)                              [4, 5, 6]
(%i4) express(u~v);
(%o4)                            [- 3, 6, - 3]
 
octave:1> u = [1 2 3]
u =

   1   2   3

octave:2> v = [4 5 6]
v =

   4   5   6

octave:3> cross(u, v)
ans =

  -3   6  -3


In [1]: import numpy as np

In [2]: u = np.array([1, 2, 3])

In [3]: v = np.array([4, 5, 6])

In [4]: np.cross(u, v)
Out[4]: array([-3,  6, -3])


 > u <- c(1, 2, 3)
 > v <- c(4, 5, 6)
 > crossprod(u, v)
      [,1]
 [1,]   32
 
$ g++ vec_prog3.cpp; ./a.out
-3.000000 6.000000 -3.000000

$ clang vec_prog3.cpp; ./a.out
-3.000000 6.000000 -3.000000

ベクトルのノルムと正規化（単位ベクトルを求める）

C言語 GLM Maxima Octave Python R言語 
#include <stdio.h>
#include <glm/glm.hpp>

int main(int argc, char *argv[]) {
	using namespace glm;
	vec3 u = vec3(1.f, 2.f, 3.f);
	{
		float ans = length(u);
		printf("%f\n", ans);
	}
	{
		vec3 ans = normalize(u);
		printf("%f %f %f\n", ans[0], ans[1], ans[2]);
	}
	return 0;
}
 
(%i1) u: [1, 2, 3];
(%o1)                              [1, 2, 3]
(%i2) sqrt(u.u), numer;
(%o2)                          3.741657386773941
(%i3) u/sqrt(u.u), numer;
(%o3)    [0.2672612419124244, 0.5345224838248488, 0.8017837257372732]

octave:1> u = [1 2 3]
u =

   1   2   3

octave:2> norm(u)
ans =  3.7417
octave:3> u/norm(u)
ans =

   0.26726   0.53452   0.80178


In [1]: import numpy as np

In [2]: u = np.array([1, 2, 3])

In [3]: np.linalg.norm(u)
Out[3]: 3.7416573867739413

In [4]: u/np.linalg.norm(u)
Out[4]: array([ 0.26726124,  0.53452248,  0.80178373])


 > u <- c(1, 2, 3)
 > sqrt(u%*%u)
          [,1]
 [1,] 3.741657
 > u/sqrt(u%*%u)
 [1] 0.2672612 0.5345225 0.8017837

$ g++ -lm vec_prog4.cpp && ./a.out
3.741657
0.267261 0.534522 0.801784

$ clang -lm vec_prog4.cpp && ./a.out
3.741657
0.267261 0.534522 0.801784

norm() 関数（や length() 関数）が無い場合は、\[ \begin{align} |\boldsymbol{u}|^2 &= \boldsymbol u\cdot \boldsymbol u \\ |\boldsymbol u| &= \sqrt{\boldsymbol u\cdot \boldsymbol u} \end{align} \]として求める

normalize() 関数が無い場合は、\( \frac{\boldsymbol u}{|\boldsymbol u|} \)とする\[ \hat{\boldsymbol u} = \frac{\boldsymbol u}{|\boldsymbol u|} \]（流儀によって、単位ベクトルを「^」（ハット）を付けて表す）

２つのベクトルのなす角を求める

C言語 GLM Maxima Octave Python R言語 
#include <stdio.h>
#include <math.h>
#include <glm/glm.hpp>

int main(int argc, char *argv[]) {
	using namespace glm;
	vec3 u = vec3(1.f, 2.f, 3.f);
	vec3 v = vec3(4.f, 5.f, 6.f);

	float ans = acosf(dot(u, v)/(length(u)*length(v)));
	printf("%f\n", ans); // radian
	printf("%f\n", ans*(180/M_PI)); // degree

	return 0;
}
 
(%i1) u: [1, 2, 3];
(%o1)                              [1, 2, 3]
(%i2) v: [4, 5, 6];
(%o2)                              [4, 5, 6]
(%i3) float(acos(u.v/(sqrt(u.u)*sqrt(v.v)))); /* radian */
(%o3)                         0.2257261285527337
(%i4) float(acos(u.v/(sqrt(u.u)*sqrt(v.v))) * (180/%pi)); /* degree */
(%o4)                          12.93315449189911

octave:1> u = [1, 2, 3]
u =

   1   2   3

octave:2> v = [4, 5, 6]
v =

   4   5   6

octave:3> acos(dot(u, v)/(norm(u)*norm(v))) % radian
ans =  0.22573
octave:4> acos(dot(u, v)/(norm(u)*norm(v))) * (180/pi) % degree
ans =  12.933

In [1]: import numpy as np

In [2]: u = np.array([1, 2, 3])

In [3]: v = np.array([4, 5, 6])

In [4]: np.arccos(u.dot(v)/(np.linalg.norm(u)*np.linalg.norm(v))) # radian
Out[4]: 0.22572612855273419

In [5]: np.arccos(u.dot(v)/(np.linalg.norm(u)*np.linalg.norm(v))) * (180/np.pi) # degree
Out[5]: 12.933154491899135


 > u <- c(1, 2, 3)
 > v <- c(4, 5, 6)
 > acos(u%*%v/(sqrt(u%*%u)*sqrt(v%*%v))) # radian
           [,1]
 [1,] 0.2257261
 > acos(u%*%v/(sqrt(u%*%u)*sqrt(v%*%v))) * (180/pi) # degree
          [,1]
 [1,] 12.93315

$ g++ -lm vec_prog5.cpp && ./a.out
0.225726
12.933170

$ clang -lm vec_prog5.cpp && ./a.out
0.225726
12.933170

内積の定義より\[ \begin{align} \boldsymbol u\cdot\boldsymbol v &= |\boldsymbol u||\boldsymbol v|\cos\theta \\ \cos\theta &= \frac{\boldsymbol u\cdot\boldsymbol v}{|\boldsymbol u||\boldsymbol v|} \\ \theta &= \arccos\frac{\boldsymbol u\cdot\boldsymbol v}{|\boldsymbol u||\boldsymbol v|} \end{align} \]

ゼロベクトル(zero vector)

ゼロベクトル \( \boldsymbol 0 \) の生成

C 言語 GLM Fortran Maxima Octave Python R 言語 
#include <stdio.h>
#include <glm/glm.hpp>

int main(int argc, char *argv[]) {
	using namespace glm;
	vec3 zero = vec3(0.f);
	printf("%f %f %f\n", zero[0], zero[1], zero[2]);

	return 0;
}
 
PROGRAM vec_prog6
  IMPLICIT NONE
  REAL(8), DIMENSION(3) :: zero = 0.D0

  WRITE(*,*) zero
END PROGRAM vec_prog6
 
(%i1) [0, 0, 0];
(%o1)                              [0, 0, 0]

octave:1> zeros(1, 3)
ans =

   0   0   0


In [1]: import numpy as np

In [2]: np.zeros(3)
Out[2]: array([ 0.,  0.,  0.])


 > rep(0, 3)
 [1] 0 0 0
 > numeric(3)
 [1] 0 0 0
 >

$ g++ vec_prog6.cpp && ./a.out
0.000000 0.000000 0.000000

$ clang vec_prog6.cpp && ./a.out
0.000000 0.000000 0.000000
 
$ gfortran vec_prog6.f90 && ./a.out
   0.0000000000000000        0.0000000000000000        0.0000000000000000

行列

行列の演算

行列とベクトルの積

C 言語 Fortran Maxima Octave Python NumPy R 言語
OpenGL Mathematics (GLM) BLAS 
#include <stdio.h>
#include <glm/vec3.hpp>
#include <glm/mat3x3.hpp>

int main(int argc, char *argv[]) {
	using namespace glm;

	mat3x3 A = mat3x3(vec3(.1f, .4f, .7f), vec3(.2f, .5f, .8f), vec3(.3f, .6f, .9f));
	printf("%f %f %f\n", A[0][0], A[1][0], A[2][0]);
	printf("%f %f %f\n", A[0][1], A[1][1], A[2][1]);
	printf("%f %f %f\n", A[0][2], A[1][2], A[2][2]);
	printf("\n");
	vec3 u = vec3(1.f, 2.f, 3.f);
	printf("%f\n%f\n%f\n", u[0], u[1], u[2]);
	printf("\n");

	vec3 ans = A*u;
	printf("%f\n%f\n%f\n", ans[0], ans[1], ans[2]);

	return 0;
}
 
#include <stdio.h>
#include <string.h>
#include <cblas.h>

int main(int argc, char *argv[]) {
	double alpha = 1.;
#define M 3
#define N 3
	double A[M][N] = {{.1, .4, .7}, {.2, .5, .8 }, { .3, .6, .9 }};
	double x[N] = { 1., 2., 3. };
	double beta = 0.;
	double y[N];
	memset(y, 0x00, sizeof(y));
#if 1 // show
	printf("%f %f %f\n", A[0][0], A[1][0], A[2][0]);
	printf("%f %f %f\n", A[0][1], A[1][1], A[2][1]);
	printf("%f %f %f\n", A[0][2], A[1][2], A[2][2]);
	printf("\n");
	{
		const double *u = x;
		printf("%f\n%f\n%f\n", u[0], u[1], u[2]);
		printf("\n");
	}
#endif

	cblas_dgemv(CblasColMajor, CblasNoTrans, M, N,
		    alpha, (double *)A, M, x, 1, beta, y, 1);

#if 1 // show
	{
		const double *ans = y;
		printf("%f\n%f\n%f\n", ans[0], ans[1], ans[2]);
	}
#endif

	return 0;
}
 
PROGRAM mat_prog1
  IMPLICIT NONE
  REAL(8), DIMENSION(3,3) :: A = RESHAPE((/.1D0, .4D0, .7D0, .2D0, .5D0, .8D0, .3D0, .6D0, .9D0/), (/3, 3/))
  REAL(8), DIMENSION(3) :: u = (/ 1., 2., 3. /)
  REAL(8), DIMENSION(3) :: ans
  INTEGER :: i, j

  DO i = 1, 3
     WRITE(*,*) (A(i,j), j = 1, 3)
  END DO
  WRITE(*,*) ''

  WRITE(*,*) u(1)
  WRITE(*,*) u(2)
  WRITE(*,*) u(3)
  WRITE(*,*) ''

  ans = MATMUL(A, u)
  WRITE(*,*) ans(1)
  WRITE(*,*) ans(2)
  WRITE(*,*) ans(3)
END PROGRAM mat_prog1
 
(%i1) A: matrix([.1, .2, .3], [.4, .5, .6], [.7, .8, .9]);
                               [ 0.1  0.2  0.3 ]
                               [               ]
(%o1)                          [ 0.4  0.5  0.6 ]
                               [               ]
                               [ 0.7  0.8  0.9 ]
(%i2) u: [1, 2, 3];
(%o2)                              [1, 2, 3]
(%i3) A.u;
                             [        1.4        ]
                             [                   ]
(%o3)                        [ 3.199999999999999 ]
                             [                   ]
                             [        5.0        ]

octave:1> A = [.1 .2 .3; .4 .5 .6; .7 .8 .9]
A =

   0.10000   0.20000   0.30000
   0.40000   0.50000   0.60000
   0.70000   0.80000   0.90000

octave:2> u = [1; 2; 3]
u =

   1
   2
   3

octave:3> A*u
ans =

   1.4000
   3.2000
   5.0000


In [1]: import numpy as np

In [2]: A = np.array([[.1, .2, .3], [.4, .5, .6], [.7, .8, .9]])

In [3]: u = np.array([[1], [2], [3]])

In [4]: A.dot(u)
Out[4]:
array([[ 1.4],
       [ 3.2],
       [ 5. ]])

In [5]: A@u
Out[5]:
array([[ 1.4],
       [ 3.2],
       [ 5. ]])


 > (A <- matrix(c(.1, .4, .7, .2, .5, .8, .3, .6, .9), ncol=3))
      [,1] [,2] [,3]
 [1,]  0.1  0.2  0.3
 [2,]  0.4  0.5  0.6
 [3,]  0.7  0.8  0.9
 > (u <- c(1, 2, 3))
 [1] 1 2 3
 > A%*%u
      [,1]
 [1,]  1.4
 [2,]  3.2
 [3,]  5.0

$ g++ mat_prog1.cpp && ./a.out
0.100000 0.200000 0.300000
0.400000 0.500000 0.600000
0.700000 0.800000 0.900000

1.000000
2.000000
3.000000

1.400000
3.200000
5.000000

$ clang mat_prog1.cpp && ./a.out
0.100000 0.200000 0.300000
0.400000 0.500000 0.600000
0.700000 0.800000 0.900000

1.000000
2.000000
3.000000

1.400000
3.200000
5.000000

$ g++ la_prog3.cpp -lblas && ./a.out
0.100000 0.200000 0.300000
0.400000 0.500000 0.600000
0.700000 0.800000 0.900000

1.000000
2.000000
3.000000

1.400000
3.200000
5.000000
 
$ gfortran mat_prog1.f90 && ./a.out
  0.10000000000000001       0.20000000000000001       0.29999999999999999     
  0.40000000000000002       0.50000000000000000       0.59999999999999998     
  0.69999999999999996       0.80000000000000004       0.90000000000000002     
 
   1.0000000000000000     
   2.0000000000000000     
   3.0000000000000000     
 
   1.3999999999999999     
   3.1999999999999997     
   5.0000000000000000

行列と行列の積

C 言語 Fortran Maxima Octave Python R 言語
GLM BLAS 
#include <stdio.h>
#include <glm/vec3.hpp>
#include <glm/mat3x3.hpp>

int main(int argc, char *argv[]) {
	using namespace glm;

	mat3x3 A = mat3x3(vec3(1.f, 4.f, -3.f), vec3(2.f, -5.f, 2.f), vec3(1.f, 6.f, 4.f));
	printf("%f %f %f\n", A[0][0], A[1][0], A[2][0]);
	printf("%f %f %f\n", A[0][1], A[1][1], A[2][1]);
	printf("%f %f %f\n", A[0][2], A[1][2], A[2][2]);
	printf("\n");
	mat3x3 B = mat3x3(vec3(1.f, 0.f, 1.f), vec3(2.f, 1.f, 1.f), vec3(-1.f, 0.f, 0.f));
	printf("%f %f %f\n", B[0][0], B[1][0], B[2][0]);
	printf("%f %f %f\n", B[0][1], B[1][1], B[2][1]);
	printf("%f %f %f\n", B[0][2], B[1][2], B[2][2]);
	printf("\n");

	mat3x3 ans = A*B;
	printf("%f %f %f\n", ans[0][0], ans[1][0], ans[2][0]);
	printf("%f %f %f\n", ans[0][1], ans[1][1], ans[2][1]);
	printf("%f %f %f\n", ans[0][2], ans[1][2], ans[2][2]);

	return 0;
}
 
#include <stdio.h>
#include <string.h>
#include <cblas.h>

int main(int argc, char *argv[]) {
	double alpha = 1.;
#define M 3
#define N 3
#define K 3
	double A[M][K] = {{1., 4., -3.}, {2., -5., 2.}, {1., 6., 4.}};
	double B[K][N] = {{1., 0., 1.}, {2., 1., 1.}, {-1., 0., 0.}};
	double beta = 0.;
	double C[M][N];
	memset(C, 0x00, sizeof(C));
#if 1 // show mat
	printf("%f %f %f\n", A[0][0], A[1][0], A[2][0]);
	printf("%f %f %f\n", A[0][1], A[1][1], A[2][1]);
	printf("%f %f %f\n", A[0][2], A[1][2], A[2][2]);
	printf("\n");
	printf("%f %f %f\n", B[0][0], B[1][0], B[2][0]);
	printf("%f %f %f\n", B[0][1], B[1][1], B[2][1]);
	printf("%f %f %f\n", B[0][2], B[1][2], B[2][2]);
	printf("\n");
#endif

	cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, K,
		    alpha, (double *)A, M, (double *)B, K, beta, (double *)C, M);
#undef M
#undef K

#if 1 // show
	{
#if 1
		const double (*ans)[N] = C;
		printf("%f %f %f\n", ans[0][0], ans[1][0], ans[2][0]);
		printf("%f %f %f\n", ans[0][1], ans[1][1], ans[2][1]);
		printf("%f %f %f\n", ans[0][2], ans[1][2], ans[2][2]);
#else
		const double *ans = (double *)C;
		printf("%f %f %f\n", ans[0], ans[3], ans[6]);
		printf("%f %f %f\n", ans[1], ans[4], ans[7]);
		printf("%f %f %f\n", ans[2], ans[5], ans[8]);
#endif
	}
#endif

#undef N

	return 0;
}
 
PROGRAM mat_prog2
  IMPLICIT NONE
  REAL(8), DIMENSION(3,3) :: A = RESHAPE((/1.D0, 4.D0, -3.D0, 2.D0, -5.D0, 2.D0, 1.D0, 6.D0, 4.D0/), (/3, 3/))
  REAL(8), DIMENSION(3,3) :: B = RESHAPE((/1.D0, 0.D0, 1.D0, 2.D0, 1.D0, 1.D0, -1.D0, 0.D0, 0.D0/), (/3, 3/))
  REAL(8), DIMENSION(3,3) :: ans
  INTEGER :: i, j

  DO i = 1, 3
     WRITE(*,*) (A(i,j), j = 1, 3)
  END DO
  WRITE(*,*) ''

  DO i = 1, 3
     WRITE(*,*) (B(i,j), j = 1, 3)
  END DO
  WRITE(*,*) ''

  ans = MATMUL(A, B)
  DO i = 1, 3
     WRITE(*,*) (ans(i,j), j = 1, 3)
  END DO
END PROGRAM mat_prog2
 
(%i1) A: matrix([1, 2, 1], [4, -5, 6], [-3, 2, 4]);
                                [  1    2   1 ]
                                [             ]
(%o1)                           [  4   - 5  6 ]
                                [             ]
                                [ - 3   2   4 ]
(%i2) B: matrix([1, 2, -1], [0, 1, 0], [1, 1, 0]);
                                 [ 1  2  - 1 ]
                                 [           ]
(%o2)                            [ 0  1   0  ]
                                 [           ]
                                 [ 1  1   0  ]
(%i3) A.B;
                                [ 2   5  - 1 ]
                                [            ]
(%o3)                           [ 10  9  - 4 ]
                                [            ]
                                [ 1   0   3  ]

octave:1> A = [1 2 1; 4 -5 6; -3 2 4]
A =

   1   2   1
   4  -5   6
  -3   2   4

octave:2> B = [1 2 -1; 0 1 0; 1 1 0]
B =

   1   2  -1
   0   1   0
   1   1   0

octave:3> A*B
ans =

    2    5   -1
   10    9   -4
    1    0    3


In [1]: import numpy as np

In [2]: A = np.array([[1, 2, 1], [4, -5, 6], [-3, 2, 4]])

In [3]: B = np.array([[1, 2, -1], [0, 1, 0], [1, 1, 0]])

In [4]: A.dot(B)
Out[4]:
array([[ 2,  5, -1],
       [10,  9, -4],
       [ 1,  0,  3]])

In [5]: A@B
Out[5]:
array([[ 2,  5, -1],
       [10,  9, -4],
       [ 1,  0,  3]])


 > (A <- matrix(c(1, 4, -3, 2, -5, 2, 1, 6, 4), ncol=3))
      [,1] [,2] [,3]
 [1,]    1    2    1
 [2,]    4   -5    6
 [3,]   -3    2    4
 > (B <- matrix(c(1, 0, 1, 2, 1, 1, -1, 0, 0), ncol=3))
      [,1] [,2] [,3]
 [1,]    1    2   -1
 [2,]    0    1    0
 [3,]    1    1    0
 > A%*%B
      [,1] [,2] [,3]
 [1,]    2    5   -1
 [2,]   10    9   -4
 [3,]    1    0    3

$ g++ mat_prog2.cpp && ./a.out
1.000000 2.000000 1.000000
4.000000 -5.000000 6.000000
-3.000000 2.000000 4.000000

1.000000 2.000000 -1.000000
0.000000 1.000000 0.000000
1.000000 1.000000 0.000000

2.000000 5.000000 -1.000000
10.000000 9.000000 -4.000000
1.000000 0.000000 3.000000

$ clang mat_prog2.cpp && ./a.out
1.000000 2.000000 1.000000
4.000000 -5.000000 6.000000
-3.000000 2.000000 4.000000

1.000000 2.000000 -1.000000
0.000000 1.000000 0.000000
1.000000 1.000000 0.000000

2.000000 5.000000 -1.000000
10.000000 9.000000 -4.000000
1.000000 0.000000 3.000000

$ g++ la_prog4.cpp -lblas && ./a.out
1.000000 2.000000 1.000000
4.000000 -5.000000 6.000000
-3.000000 2.000000 4.000000

1.000000 2.000000 -1.000000
0.000000 1.000000 0.000000
1.000000 1.000000 0.000000

2.000000 5.000000 -1.000000
10.000000 9.000000 -4.000000
1.000000 0.000000 3.000000
 
$ gfortran mat_prog2.f90 && ./a.out
   1.0000000000000000        2.0000000000000000        1.0000000000000000     
   4.0000000000000000       -5.0000000000000000        6.0000000000000000     
  -3.0000000000000000        2.0000000000000000        4.0000000000000000     
 
   1.0000000000000000        2.0000000000000000       -1.0000000000000000     
   0.0000000000000000        1.0000000000000000        0.0000000000000000     
   1.0000000000000000        1.0000000000000000        0.0000000000000000     
 
   2.0000000000000000        5.0000000000000000       -1.0000000000000000     
   10.000000000000000        9.0000000000000000       -4.0000000000000000     
   1.0000000000000000        0.0000000000000000        3.0000000000000000

（Pythonでは、Python 3.5から「@」演算子が提供されている）

転置行列

C言語 GLM Fortran Maxima Octave Python R言語 
#include <stdio.h>
#include <glm/mat3x3.hpp>
#include <glm/matrix.hpp>

int main(int argc, char *argv[]) {
	using namespace glm;

	mat3x3 A = mat3x3(vec3(1.f, 4.f, 7.f), vec3(2.f, 5.f, 8.f), vec3(3.f, 6.f, 9.f));
	printf("%f %f %f\n", A[0][0], A[1][0], A[2][0]);
	printf("%f %f %f\n", A[0][1], A[1][1], A[2][1]);
	printf("%f %f %f\n", A[0][2], A[1][2], A[2][2]);
	printf("\n");
	
	mat3x3 ans = transpose(A);
	printf("%f %f %f\n", ans[0][0], ans[1][0], ans[2][0]);
	printf("%f %f %f\n", ans[0][1], ans[1][1], ans[2][1]);
	printf("%f %f %f\n", ans[0][2], ans[1][2], ans[2][2]);

	return 0;
}
 
PROGRAM mat_prog3
  IMPLICIT NONE
  REAL(8), DIMENSION(3,3) :: A = RESHAPE((/1.D0, 4.D0, 7.D0, 2.D0, 5.D0, 8.D0, 3.D0, 6.D0, 9.D0/), (/3, 3/))
  REAL(8), DIMENSION(3,3) :: ans
  INTEGER :: i, j

  DO i = 1, 3
     WRITE(*,*) (A(i,j), j = 1, 3)
  END DO
  WRITE(*,*) ''

  ans = TRANSPOSE(A)
  DO i = 1, 3
     WRITE(*,*) (ans(i,j), j = 1, 3)
  END DO
END PROGRAM mat_prog3
 
(%i1) A: matrix([1, 2, 3], [4, 5, 6], [7, 8, 9]);
                                  [ 1  2  3 ]
                                  [         ]
(%o1)                             [ 4  5  6 ]
                                  [         ]
                                  [ 7  8  9 ]
(%i2) transpose(A);
                                  [ 1  4  7 ]
                                  [         ]
(%o2)                             [ 2  5  8 ]
                                  [         ]
                                  [ 3  6  9 ]

octave:1> A = [1 2 3; 4 5 6; 7 8 9]
A =

   1   2   3
   4   5   6
   7   8   9

octave:2> A'
ans =

   1   4   7
   2   5   8
   3   6   9


In [1]: import numpy as np

In [2]: A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

In [3]: A.T
Out[3]:
array([[1, 4, 7],
       [2, 5, 8],
       [3, 6, 9]])


 > (A <- matrix(c(1, 4, 7, 2, 5, 8, 3, 6, 9), ncol=3))
      [,1] [,2] [,3]
 [1,]    1    2    3
 [2,]    4    5    6
 [3,]    7    8    9
 > t(A)
      [,1] [,2] [,3]
 [1,]    1    4    7
 [2,]    2    5    8
 [3,]    3    6    9

$ g++ mat_prog3.cpp && ./a.out
1.000000 2.000000 3.000000
4.000000 5.000000 6.000000
7.000000 8.000000 9.000000

1.000000 4.000000 7.000000
2.000000 5.000000 8.000000
3.000000 6.000000 9.000000

$ clang mat_prog3.cpp && ./a.out
1.000000 2.000000 3.000000
4.000000 5.000000 6.000000
7.000000 8.000000 9.000000

1.000000 4.000000 7.000000
2.000000 5.000000 8.000000
3.000000 6.000000 9.000000

$ gfortran mat_prog3.f90 && ./a.out
   1.0000000000000000        2.0000000000000000        3.0000000000000000     
   4.0000000000000000        5.0000000000000000        6.0000000000000000     
   7.0000000000000000        8.0000000000000000        9.0000000000000000     
 
   1.0000000000000000        4.0000000000000000        7.0000000000000000     
   2.0000000000000000        5.0000000000000000        8.0000000000000000     
   3.0000000000000000        6.0000000000000000        9.0000000000000000

（Octave では、「'」は共役転置行列。ただの転置行列は「.'」）

逆行列

C言語 GLM Maxima Octave Python R言語 
#include <stdio.h>
#include <glm/mat3x3.hpp>
#include <glm/matrix.hpp>

int main(int argc, char *argv[]) {
	using namespace glm;

	mat3x3 A = mat3x3(vec3(-1.f, 4.f, 7.f), vec3(2.f, -5.f, 8.f), vec3(3.f, 6.f, -9.f));
	printf("%f %f %f\n", A[0][0], A[1][0], A[2][0]);
	printf("%f %f %f\n", A[0][1], A[1][1], A[2][1]);
	printf("%f %f %f\n", A[0][2], A[1][2], A[2][2]);
	printf("\n");

	mat3x3 ans = inverse(A);
	printf("%f %f %f\n", ans[0][0], ans[1][0], ans[2][0]);
	printf("%f %f %f\n", ans[0][1], ans[1][1], ans[2][1]);
	printf("%f %f %f\n", ans[0][2], ans[1][2], ans[2][2]);

	return 0;
}
 
(%i1) A: matrix([-1, 2, 3], [4, -5, 6], [7, 8, -9]);
                               [ - 1   2    3  ]
                               [               ]
(%o1)                          [  4   - 5   6  ]
                               [               ]
                               [  7    8   - 9 ]
(%i2) invert(A), numer;
      [ - 0.008333333333333334   0.1166666666666666            0.075          ]
      [                                                                       ]
(%o2) [   0.2166666666666666    - 0.03333333333333333           0.05          ]
      [                                                                       ]
      [   0.1861111111111111     0.06111111111111111   - 0.008333333333333334 ]
(%i3) A^^-1, numer;
      [ - 0.008333333333333334   0.1166666666666666            0.075          ]
      [                                                                       ]
(%o3) [   0.2166666666666666    - 0.03333333333333333           0.05          ]
      [                                                                       ]
      [   0.1861111111111111     0.06111111111111111   - 0.008333333333333334 ]

octave:1> A = [-1 2 3; 4 -5 6; 7 8 -9]
A =

  -1   2   3
   4  -5   6
   7   8  -9

octave:2> inv(A)
ans =

  -0.0083333   0.1166667   0.0750000
   0.2166667  -0.0333333   0.0500000
   0.1861111   0.0611111  -0.0083333

octave:3> A^-1
ans =

  -0.0083333   0.1166667   0.0750000
   0.2166667  -0.0333333   0.0500000
   0.1861111   0.0611111  -0.0083333

octave:4> A**-1
ans =

  -0.0083333   0.1166667   0.0750000
   0.2166667  -0.0333333   0.0500000
   0.1861111   0.0611111  -0.0083333


In [1]: import numpy as np

In [2]: A = np.array([[-1, 2, 3], [4, -5, 6], [7, 8, -9]])

In [3]: np.linalg.inv(A)
Out[3]:
array([[-0.00833333,  0.11666667,  0.075     ],
       [ 0.21666667, -0.03333333,  0.05      ],
       [ 0.18611111,  0.06111111, -0.00833333]])


 > (A <- matrix(c(-1, 4, 7, 2, -5, 8, 3, 6, -9), ncol=3))
      [,1] [,2] [,3]
 [1,]   -1    2    3
 [2,]    4   -5    6
 [3,]    7    8   -9
 > solve(A)
              [,1]        [,2]         [,3]
 [1,] -0.008333333  0.11666667  0.075000000
 [2,]  0.216666667 -0.03333333  0.050000000
 [3,]  0.186111111  0.06111111 -0.008333333

$ g++ mat_prog4.cpp && ./a.out
-1.000000 2.000000 3.000000
4.000000 -5.000000 6.000000
7.000000 8.000000 -9.000000

-0.008333 0.116667 0.075000
0.216667 -0.033333 0.050000
0.186111 0.061111 -0.008333

$ clang mat_prog4.cpp && ./a.out
-1.000000 2.000000 3.000000
4.000000 -5.000000 6.000000
7.000000 8.000000 -9.000000

-0.008333 0.116667 0.075000
0.216667 -0.033333 0.050000
0.186111 0.061111 -0.008333

単位行列

C 言語 GLM Fortran Maxima Octave Python R 言語 
#include <stdio.h>
#include <glm/mat3x3.hpp>

int main(int argc, char *argv[]) {
	using namespace glm;

	mat3x3 I = mat3x3(1.f);
	printf("%f %f %f\n", I[0][0], I[1][0], I[2][0]);
	printf("%f %f %f\n", I[0][1], I[1][1], I[2][1]);
	printf("%f %f %f\n", I[0][2], I[1][2], I[2][2]);

	return 0;
}
 
PROGRAM mat_prog5
  IMPLICIT NONE
  INTEGER :: i, j
  REAL(8), DIMENSION(3,3) :: D = RESHAPE((/1.D0, (0.D0, i = 1, 3), 1.D0, (0.D0, i = 1, 3), 1.D0/), (/3, 3/))

  DO i = 1, 3
     WRITE(*,*) (D(i,j), j = 1, 3)
  END DO
END PROGRAM mat_prog5
 
(%i1) ident(3);
                                  [ 1  0  0 ]
                                  [         ]
(%o1)                             [ 0  1  0 ]
                                  [         ]
                                  [ 0  0  1 ]

octave:1> eye(3)
ans =

Diagonal Matrix

   1   0   0
   0   1   0
   0   0   1


In [1]: import numpy as np

In [2]: np.eye(3)
Out[2]:
array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])


 > diag(rep(1, 3))
      [,1] [,2] [,3]
 [1,]    1    0    0
 [2,]    0    1    0
 [3,]    0    0    1

$ g++ mat_prog5.cpp && ./a.out
1.000000 0.000000 0.000000
0.000000 1.000000 0.000000
0.000000 0.000000 1.000000

$ clang mat_prog5.cpp && ./a.out
1.000000 0.000000 0.000000
0.000000 1.000000 0.000000
0.000000 0.000000 1.000000

$ gfortran mat_prog5.f90 && ./a.out
   1.0000000000000000        0.0000000000000000        0.0000000000000000     
   0.0000000000000000        1.0000000000000000        0.0000000000000000     
   0.0000000000000000        0.0000000000000000        1.0000000000000000

スカラーと行列の積

Fortran Maxima Octave Python R 言語
PROGRAM mat_prog6
  IMPLICIT NONE
  INTEGER :: i, j
  REAL(8), DIMENSION(3,3) :: E = RESHAPE((/1.D0, (0.D0, i = 1, 3), 1.D0, (0.D0, i = 1, 3), 1.D0/), (/3, 3/))
  REAL(8), DIMENSION(3,3) :: ans

  ans = 2.D0*E
  DO i = 1, 3
     WRITE(*,*) (ans(i,j), j = 1, 3)
  END DO
END PROGRAM mat_prog6
 
(%i1) 2*ident(3);
                                  [ 2  0  0 ]
                                  [         ]
(%o1)                             [ 0  2  0 ]
                                  [         ]
                                  [ 0  0  2 ]

octave:1> 2*eye(3)
ans =

Diagonal Matrix

   2   0   0
   0   2   0
   0   0   2


In [1]: import numpy as np

In [2]: 2*np.eye(3)
Out[2]:
array([[ 2.,  0.,  0.],
       [ 0.,  2.,  0.],
       [ 0.,  0.,  2.]])


 > 2*diag(3)
      [,1] [,2] [,3]
 [1,]    2    0    0
 [2,]    0    2    0
 [3,]    0    0    2

$ gfortran mat_prog6.f90 && ./a.out
   2.0000000000000000        0.0000000000000000        0.0000000000000000     
   0.0000000000000000        2.0000000000000000        0.0000000000000000     
   0.0000000000000000        0.0000000000000000        2.0000000000000000

行列のトレースを求める

Maxima Octave Python R言語
(%i1) load("nchrpl");
(%o1)          /usr/share/maxima/5.37.2/share/matrix/nchrpl.mac
(%i2) A: matrix([1, 0, 0], [0, 2, 0], [0, 0, 3]);
                                  [ 1  0  0 ]
                                  [         ]
(%o2)                             [ 0  2  0 ]
                                  [         ]
                                  [ 0  0  3 ]
(%i3) mattrace(A);
(%o3)                                  6

octave:1> A = diag([1 2 3])
A =

Diagonal Matrix

   1   0   0
   0   2   0
   0   0   3

octave:2> trace(A)
ans =  6

In [1]: import numpy as np

In [2]: A = np.diag([1, 2, 3])

In [3]: A.trace()
Out[3]: 6


 > (A <- diag(c(1, 2, 3)))
      [,1] [,2] [,3]
 [1,]    1    0    0
 [2,]    0    2    0
 [3,]    0    0    3
 > sum(diag(A))
 [1] 6

行列のランクを求める

Maxima Octave Python R 言語
(%i1) A: matrix([-1, 2, 3], [4, -5, 6], [7, 8, -9]);
                               [ - 1   2    3  ]
                               [               ]
(%o1)                          [  4   - 5   6  ]
                               [               ]
                               [  7    8   - 9 ]
(%i2) rank(A);
(%o2)                                  3

octave:1> A = [-1 2 3; 4 -5 6; 7 8 -9]
A =

  -1   2   3
   4  -5   6
   7   8  -9

octave:2> rank(A)
ans =  3

In [1]: import numpy as np

In [2]: A = np.array([[-1, 2, 3], [4, -5, 6], [7, 8, -9]])

In [3]: np.linalg.matrix_rank(A)
Out[3]: 3


 > (A <- matrix(c(-1, 4, 7, 2, -5, 8, 3, 6, -9), ncol=3))
      [,1] [,2] [,3]
 [1,]   -1    2    3
 [2,]    4   -5    6
 [3,]    7    8   -9
 > qr(A)$rank
 [1] 3

行列式の値を求める

C 言語 GLM Maxima Octave Python R 言語
#include <stdio.h>
#include <glm/vec3.hpp>
#include <glm/mat3x3.hpp>
#include <glm/matrix.hpp>

int main(int argc, char *argv[]) {
	using namespace glm;

	mat3x3 A(vec3(-1.f, 4.f, 7.f), vec3(2.f, -5.f, 8.f), vec3(3.f, 6.f, -9.f));
	printf("%f %f %f\n", A[0][0], A[1][0], A[2][0]);
	printf("%f %f %f\n", A[0][1], A[1][1], A[2][1]);
	printf("%f %f %f\n", A[0][2], A[1][2], A[2][2]);
	printf("\n");

	float ans = determinant(A);
	printf("%f\n", ans);
}

$ g++ mat_prog6.cpp && ./a.out
-1.000000 2.000000 3.000000
4.000000 -5.000000 6.000000
7.000000 8.000000 -9.000000

360.000000

(%i1) A: matrix([-1, 2, 3], [4, -5, 6], [7, 8, -9]);
                               [ - 1   2    3  ]
                               [               ]
(%o1)                          [  4   - 5   6  ]
                               [               ]
                               [  7    8   - 9 ]
(%i2) determinant(A);
(%o2)                                 360

octave:1> A = [-1 2 3; 4 -5 6; 7 8 -9]
A =

  -1   2   3
   4  -5   6
   7   8  -9

octave:2> det(A)
ans =  360.00

In [1]: import numpy as np

In [2]: A = np.array([[-1, 2, 3], [4, -5, 6], [7, 8, -9]])

In [3]: np.linalg.det(A)
Out[3]: 360.00000000000006


 > (A <- matrix(c(-1, 4, 7, 2, -5, 8, 3, 6, -9), ncol=3))
      [,1] [,2] [,3]
 [1,]   -1    2    3
 [2,]    4   -5    6
 [3,]    7    8   -9
 > det(A)
 [1] 360

連立一次方程式の解法

係数行列 \( A \) と定数ベクトル \( \boldsymbol b \) を用いて解ベクトル \( \boldsymbol x \) を求める\[ \begin{align} A\boldsymbol x &= \boldsymbol b \\ \boldsymbol x &= A^{-1}\boldsymbol b \end{align} \]

C言語 LAPACK Octave Python R言語 
#include <stdio.h>
#include <string.h>
#include <lapacke.h>

int main(int argc, char *argv[]) {
#define N 3
	double A[N][N] = {{-1., 2., 3.}, {4., -5., 6.}, {7., 8., -9.}};
	double b[N] = {1., 1., 1.};
#if 1 // show A, b
	printf("%f %f %f\n", A[0][0], A[0][1], A[0][2]);
	printf("%f %f %f\n", A[1][0], A[1][1], A[1][2]);
	printf("%f %f %f\n", A[2][0], A[2][1], A[2][2]);
	printf("\n");
	printf("%f\n%f\n%f\n", b[0], b[1], b[2]);
	printf("\n");
#endif
	lapack_int pivot[N];
	memset(pivot, 0x00, sizeof(pivot));

	lapack_int info =
		LAPACKE_dgesv(LAPACK_ROW_MAJOR, N, 1, (double *)A, N, pivot, b, 1);
	if (info < 0)
		printf("RETURN VALUE: %d\n", info);
#if 1 // show x
	{
		const double *x = b;
		printf("%f\n%f\n%f\n", x[0], x[1], x[2]);
		printf("\n");
	}
#endif
	return 0;
}

octave:1> A = [-1 2 3; 4 -5 6; 7 8 -9]
A =

  -1   2   3
   4  -5   6
   7   8  -9

octave:2> b = [1; 1; 1]
b =

   1
   1
   1

octave:3> x = A\b
x =

   0.18333
   0.23333
   0.23889


In [1]: import numpy as np

In [2]: A = np.array([[-1, 2, 3], [4, -5, 6], [7, 8, -9]]); A
Out[2]:
array([[-1,  2,  3],
       [ 4, -5,  6],
       [ 7,  8, -9]])

In [3]: b = [[1], [1], [1]]; b
Out[3]: [[1], [1], [1]]

In [4]: x = np.linalg.solve(A, b); x
Out[4]:
array([[ 0.18333333],
       [ 0.23333333],
       [ 0.23888889]])


 > (A <- matrix(c(-1, 4, 7, 2, -5, 8, 3, 6, -9), ncol=3))
      [,1] [,2] [,3]
 [1,]   -1    2    3
 [2,]    4   -5    6
 [3,]    7    8   -9
 > (b <- rep(1, 3))
 [1] 1 1 1
 > (x <- solve(A, b))
 [1] 0.1833333 0.2333333 0.2388889

$ g++ la_prog5.cpp -llapacke && ./a.out
-1.000000 2.000000 3.000000
4.000000 -5.000000 6.000000
7.000000 8.000000 -9.000000

1.000000
1.000000
1.000000

0.183333
0.233333
0.238889

Octave では、左除算の演算子「\( A\backslash b \)」が用意されている

Octave 
octave:1> 2\1
ans =  0.50000
「\( \backslash \)」演算子で、逆行列を経由せず「ガウスの消去法」で解いている\[ \begin{align} A\boldsymbol x &= \boldsymbol b \\ \boldsymbol x &= A^{-1}\boldsymbol b \\ \boldsymbol x &= A \backslash\boldsymbol b \end{align} \]（連立方程式を解くために逆行列を使うのは計算コストの無駄であり、さらに、計算精度が悪くなる）

Maximaでは、方程式と変数をリストで渡す

Maxima 
(%i1) linsolve([-1*x + 2*y + 3*z = 1, 4*x + -5*y + 6*z = 1, 7*x + 8*y + -9*z = 1], [x, y, z]), numer;
(%o1) [x = 0.1833333333333333, y = 0.2333333333333333, z = 0.2388888888888889]

係数行列、拡大係数行列(augmented matrix)を取り出す

Maxima 
(%i1) A: coefmatrix([-1*x + 2*y + 3*z = 1, 4*x -5*y + 6*z = 1, 7*x + 8*y -9*z = 1], [x, y, z]);
                               [ - 1   2    3  ]
                               [               ]
(%o1)                          [  4   - 5   6  ]
                               [               ]
                               [  7    8   - 9 ]
(%i2) A: augcoefmatrix([-1*x + 2*y + 3*z = 1, 4*x -5*y + 6*z = 1, 7*x + 8*y -9*z = 1], [x, y, z]);
                            [ - 1   2    3   - 1 ]
                            [                    ]
(%o2)                       [  4   - 5   6   - 1 ]
                            [                    ]
                            [  7    8   - 9  - 1 ]

（Maxima の augcoefmatrix() 関数は、\( (A|\boldsymbol b) \) ではなく \( (A|{-\boldsymbol b}) \) を求める）

乱数行列

Octave Python R言語 
octave:4> rand("seed", 0)
octave:5> A = rand(4, 4)
A =

   1.000000   0.036945   0.204136   0.910082
   0.974520   0.161712   0.167307   0.111716
   0.647484   0.664645   0.654957   0.299021
   0.333086   0.084673   0.128816   0.265447

octave:6> b = rand(4, 1)
b =

   0.699824
   0.950282
   0.091811
   0.902341

octave:7> x = A\b
x =

     0.94201
  -237.55854
   262.92529
   -49.59809


In [5]: np.random.seed(0)

In [6]: A = np.random.rand(4, 4); A
Out[6]:
array([[ 0.5488135 ,  0.71518937,  0.60276338,  0.54488318],
       [ 0.4236548 ,  0.64589411,  0.43758721,  0.891773  ],
       [ 0.96366276,  0.38344152,  0.79172504,  0.52889492],
       [ 0.56804456,  0.92559664,  0.07103606,  0.0871293 ]])

In [7]: b = np.random.rand(4, 1); b
Out[7]:
array([[ 0.0202184 ],
       [ 0.83261985],
       [ 0.77815675],
       [ 0.87001215]])

In [8]: x = np.linalg.solve(A, b); x
Out[8]:
array([[ 2.80445897],
       [-0.69248685],
       [-3.21927101],
       [ 1.68258328]])


 > set.seed(0)
 > (A <- matrix(runif(4*4), ncol=4))
           [,1]      [,2]       [,3]      [,4]
 [1,] 0.8966972 0.9082078 0.66079779 0.1765568
 [2,] 0.2655087 0.2016819 0.62911404 0.6870228
 [3,] 0.3721239 0.8983897 0.06178627 0.3841037
 [4,] 0.5728534 0.9446753 0.20597457 0.7698414
 > (b <- runif(4))
 [1] 0.4976992 0.7176185 0.9919061 0.3800352
 > (x <- solve(A, b))
 [1] -4.4439116  3.0436875  2.7815924 -0.6786966

疑似逆行列

Octave Python R言語
octave:1> A = [1 2 3 4; 5 6 7 8]
A =

   1   2   3   4
   5   6   7   8

octave:2> pinv(A)
ans =

  -5.5000e-01   2.5000e-01
  -2.2500e-01   1.2500e-01
   1.0000e-01  -1.3878e-17
   4.2500e-01  -1.2500e-01


In [1]: import numpy as np

In [2]: A = np.array([[1, 2, 3, 4], [5, 6, 7, 8]]); A
Out[2]:
array([[1, 2, 3, 4],
       [5, 6, 7, 8]])

In [3]: np.linalg.pinv(A)
Out[3]:
array([[ -5.50000000e-01,   2.50000000e-01],
       [ -2.25000000e-01,   1.25000000e-01],
       [  1.00000000e-01,  -2.08166817e-17],
       [  4.25000000e-01,  -1.25000000e-01]])


 > library(MASS)
 > (A <- matrix(c(1, 5, 2, 6, 3, 7, 4, 8), ncol=4))
      [,1] [,2] [,3] [,4]
 [1,]    1    2    3    4
 [2,]    5    6    7    8
 > ginv(A)
        [,1]          [,2]
 [1,] -0.550  2.500000e-01
 [2,] -0.225  1.250000e-01
 [3,]  0.100 -2.081668e-17
 [4,]  0.425 -1.250000e-01

固有値と固有ベクトル
Maxima Octave Python R言語
(%i1) load("eigen");
(%o1)           /usr/share/maxima/5.37.2/share/matrix/eigen.mac
(%i2) A: matrix([1, 2], [3, 4]);
                                   [ 1  2 ]
(%o2)                              [      ]
                                   [ 3  4 ]
(%i3) float(uniteigenvectors(A));
(%o3) [[[- 0.3722813232690143, 5.372281323269014], [1.0, 1.0]],
[[[0.824564840132394, - 0.5657674649689922]],
[[0.4159735579192844, 0.9093767091321241]]]]

octave:1> A = [1 2; 3 4]
A =

   1   2
   3   4

octave:2> [V, D] = eig(A)
V =

  -0.82456  -0.41597
   0.56577  -0.90938

D =

Diagonal Matrix

  -0.37228         0
         0   5.37228


In [1]: import numpy as np

In [2]: A = np.array([[1, 2], [3, 4]]); A
Out[2]:
array([[1, 2],
       [3, 4]])

In [3]: val, vec = np.linalg.eig(A)

In [4]: val
Out[4]: array([-0.37228132,  5.37228132])

In [5]: vec
Out[5]:
array([[-0.82456484, -0.41597356],
       [ 0.56576746, -0.90937671]])


 > (A <- matrix(c(1, 3, 2, 4), ncol=2))
      [,1] [,2]
 [1,]    1    2
 [2,]    3    4
 > eigen(A)
 $values
 [1]  5.3722813 -0.3722813

 $vectors
            [,1]       [,2]
 [1,] -0.4159736 -0.8245648
 [2,] -0.9093767  0.5657675

特異値分解 (SVD, singular value decomposition)

\[ A = U\Sigma V^\mathrm T \]

Octave 
octave:1> rand("seed", 0)
octave:2> rand(1); % adhoc
octave:3> A = [pascal(3), rand(3, 1)]
A =

   1.00000   1.00000   1.00000   0.97452
   1.00000   2.00000   3.00000   0.64748
   1.00000   3.00000   6.00000   0.33309

octave:4> [U S V] = svd(A)
U =

  -0.20438   0.84690   0.49090
  -0.47582   0.35231  -0.80590
  -0.85547  -0.39829   0.33097

S =

Diagonal Matrix

   7.91219         0         0         0
         0   1.36110         0         0
         0         0   0.15680         0

V =

  -0.1940887   0.5884353   0.1018183  -0.7782715
  -0.4704668   0.2620311  -0.8163762   0.2086401
  -0.8549654  -0.3569951   0.3762084  -0.0074839
  -0.1001241   0.6764896   0.4261810   0.5922050

（Octave では、pascal() 関数は「Pascal 行列」を返す）

© 2020 Tadakazu Nagai