ベクトルと行列

永井 忠一 2020.5.30

ベクトル

ベクトルの演算

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

C言語 OpenGL Mathematics (GLM) Maxima Octave Python NumPy R言語 
#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;
}
 
(%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
C言語 BLAS
#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;
}
$ 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

ベクトルの内積(dot 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);

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

	return 0;
}
 
(%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
C言語 BLAS
#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;
}
$ g++ la_prog2.cpp -lblas && ./a.out
32.000000

ベクトルの外積(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 <math.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 = sqrtf(dot(u, 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() 関数が無い場合は、\[ \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)/(sqrtf(dot(u, u))*sqrtf(dot(v, 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 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;
}
 
(%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

行列

行列の演算

行列とベクトルの積

C言語 OpenGL Mathematics (GLM) Maxima Octave Python NumPy R言語 
#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;
}
 
(%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
C言語 BLAS
#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;
}
$ 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

行列と行列の積

C言語 GLM Maxima Octave Python R言語 
#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;
}
 
(%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
C言語 BLAS
#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;
}
$ 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

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

転置行列

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 = 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;
}
 
(%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

（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 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;
}
 
(%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

スカラーと行列の積

Maxima Octave Python R言語
(%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

行列のトレースを求める

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
(%i3) determinant(A);
(%o3)                                 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> rank(A)
ans =  3
octave:3> 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.matrix_rank(A)
Out[3]: 3

In [4]: np.linalg.det(A)
Out[4]: 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
 > qr(A)$rank
 [1] 3
 > 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