diff --git a/colorm.go b/colorm.go
index cd6a56810..57be2f31e 100644
--- a/colorm.go
+++ b/colorm.go
@@ -36,7 +36,7 @@ type ColorM struct {
 // Concat multiplies a color matrix with the other color matrix.
 // This is same as muptiplying the matrix other and the matrix c in this order.
 func (c *ColorM) Concat(other ColorM) {
-	c.impl.Concat(other.impl)
+	c.impl.Concat(&other.impl)
 }
 
 // Add is deprecated as of 1.5.0-alpha.
diff --git a/geom.go b/geom.go
index 75f68c0b3..fb54a16ee 100644
--- a/geom.go
+++ b/geom.go
@@ -30,13 +30,29 @@ type GeoM struct {
 
 // Element returns a value of a matrix at (i, j).
 func (g *GeoM) Element(i, j int) float64 {
-	return g.impl.UnsafeElements()[i*affine.GeoMDim+j]
+	a, b, c, d, tx, ty := g.impl.Elements()
+	switch {
+	case i == 0 && j == 0:
+		return a
+	case i == 0 && j == 1:
+		return b
+	case i == 0 && j == 2:
+		return tx
+	case i == 1 && j == 0:
+		return c
+	case i == 1 && j == 1:
+		return d
+	case i == 1 && j == 2:
+		return ty
+	default:
+		panic("ebiten: i or j is out of index")
+	}
 }
 
 // Concat multiplies a geometry matrix with the other geometry matrix.
 // This is same as muptiplying the matrix other and the matrix g in this order.
 func (g *GeoM) Concat(other GeoM) {
-	g.impl.Concat(other.impl)
+	g.impl.Concat(&other.impl)
 }
 
 // Add is deprecated as of 1.5.0-alpha.
diff --git a/internal/affine/colorm.go b/internal/affine/colorm.go
index 2519e2d0a..c14632091 100644
--- a/internal/affine/colorm.go
+++ b/internal/affine/colorm.go
@@ -82,7 +82,7 @@ func (c *ColorM) Equals(other *ColorM) bool {
 
 // Concat multiplies a color matrix with the other color matrix.
 // This is same as muptiplying the matrix other and the matrix c in this order.
-func (c *ColorM) Concat(other ColorM) {
+func (c *ColorM) Concat(other *ColorM) {
 	if c.elements == nil {
 		c.elements = colorMIdentityElements
 	}
@@ -182,8 +182,8 @@ var (
 // This conversion uses RGB to/from YCrCb conversion.
 func (c *ColorM) ChangeHSV(hueTheta float64, saturationScale float64, valueScale float64) {
 	sin, cos := math.Sincos(hueTheta)
-	c.Concat(rgbToYCbCr)
-	c.Concat(ColorM{
+	c.Concat(&rgbToYCbCr)
+	c.Concat(&ColorM{
 		elements: []float64{
 			1, 0, 0, 0, 0,
 			0, cos, -sin, 0, 0,
@@ -194,7 +194,7 @@ func (c *ColorM) ChangeHSV(hueTheta float64, saturationScale float64, valueScale
 	s := saturationScale
 	v := valueScale
 	c.Scale(v, s*v, s*v, 1)
-	c.Concat(yCbCrToRgb)
+	c.Concat(&yCbCrToRgb)
 }
 
 var monochrome ColorM
diff --git a/internal/affine/geom.go b/internal/affine/geom.go
index 79c1de415..3cfa30a04 100644
--- a/internal/affine/geom.go
+++ b/internal/affine/geom.go
@@ -21,105 +21,138 @@ import (
 // GeoMDim is a dimension of a GeoM.
 const GeoMDim = 3
 
-var (
-	geoMIdentityElements = []float64{
-		1, 0, 0,
-		0, 1, 0,
-	}
-)
-
 // A GeoM represents a matrix to transform geometry when rendering an image.
 //
 // The initial value is identity.
 type GeoM struct {
-	// When elements is empty, this matrix is identity.
-	// elements is immutable and a new array must be created when updating.
-	elements []float64
+	a      float64
+	b      float64
+	c      float64
+	d      float64
+	tx     float64
+	ty     float64
+	inited bool
 }
 
-func (g *GeoM) UnsafeElements() []float64 {
-	if g.elements == nil {
-		g.elements = geoMIdentityElements
+func (g *GeoM) Elements() (a, b, c, d, tx, ty float64) {
+	if !g.inited {
+		return 1, 0, 0, 1, 0, 0
 	}
-	return g.elements
+	return g.a, g.b, g.c, g.d, g.tx, g.ty
+}
+
+func (g *GeoM) init() {
+	g.a = 1
+	g.b = 0
+	g.c = 0
+	g.d = 1
+	g.tx = 0
+	g.ty = 0
+	g.inited = true
 }
 
 // SetElement sets an element at (i, j).
 func (g *GeoM) SetElement(i, j int, element float64) {
-	if g.elements == nil {
-		g.elements = geoMIdentityElements
+	if !g.inited {
+		g.init()
+	}
+	switch {
+	case i == 0 && j == 0:
+		g.a = element
+	case i == 0 && j == 1:
+		g.b = element
+	case i == 0 && j == 2:
+		g.tx = element
+	case i == 1 && j == 0:
+		g.c = element
+	case i == 1 && j == 1:
+		g.d = element
+	case i == 1 && j == 2:
+		g.ty = element
+	default:
+		panic("affine: i or j is out of index")
 	}
-	es := make([]float64, len(g.elements))
-	copy(es, g.elements)
-	es[i*GeoMDim+j] = element
-	g.elements = es
 }
 
 // Concat multiplies a geometry matrix with the other geometry matrix.
 // This is same as muptiplying the matrix other and the matrix g in this order.
-func (g *GeoM) Concat(other GeoM) {
-	if g.elements == nil {
-		g.elements = geoMIdentityElements
+func (g *GeoM) Concat(other *GeoM) {
+	if !g.inited {
+		g.init()
 	}
-	if other.elements == nil {
-		other.elements = geoMIdentityElements
+	if !other.inited {
+		other.init()
 	}
-	g.elements = mul(other.elements, g.elements, GeoMDim)
+	a, b, c, d, tx, ty := g.a, g.b, g.c, g.d, g.tx, g.ty
+	g.a = other.a*a + other.b*c
+	g.b = other.a*b + other.b*d
+	g.tx = other.a*tx + other.b*ty + other.tx
+	g.c = other.c*a + other.d*c
+	g.d = other.c*b + other.d*d
+	g.ty = other.c*tx + other.d*ty + other.ty
 }
 
 // Add is deprecated.
 func (g *GeoM) Add(other GeoM) {
-	if g.elements == nil {
-		g.elements = geoMIdentityElements
+	if !g.inited {
+		g.init()
 	}
-	if other.elements == nil {
-		other.elements = geoMIdentityElements
+	if !other.inited {
+		other.init()
 	}
-	g.elements = add(other.elements, g.elements, GeoMDim)
+	g.a += other.a
+	g.b += other.b
+	g.c += other.c
+	g.d += other.d
+	g.tx += other.tx
+	g.ty += other.ty
 }
 
 // Scale scales the matrix by (x, y).
 func (g *GeoM) Scale(x, y float64) {
-	if g.elements == nil {
-		g.elements = []float64{
-			x, 0, 0,
-			0, y, 0,
-		}
+	if !g.inited {
+		g.a = x
+		g.b = 0
+		g.c = 0
+		g.d = y
+		g.tx = 0
+		g.ty = 0
+		g.inited = true
 		return
 	}
-	es := make([]float64, len(g.elements))
-	copy(es, g.elements)
-	for i := 0; i < GeoMDim; i++ {
-		es[i] *= x
-		es[i+GeoMDim] *= y
-	}
-	g.elements = es
+	g.a *= x
+	g.b *= x
+	g.tx *= x
+	g.c *= y
+	g.d *= y
+	g.ty *= y
 }
 
 // Translate translates the matrix by (x, y).
 func (g *GeoM) Translate(tx, ty float64) {
-	if g.elements == nil {
-		g.elements = []float64{
-			1, 0, tx,
-			0, 1, ty,
-		}
+	if !g.inited {
+		g.a = 1
+		g.b = 0
+		g.c = 0
+		g.d = 1
+		g.tx = tx
+		g.ty = ty
+		g.inited = true
 		return
 	}
-	es := make([]float64, len(g.elements))
-	copy(es, g.elements)
-	es[2] += tx
-	es[2+GeoMDim] += ty
-	g.elements = es
+	g.tx += tx
+	g.ty += ty
 }
 
 // Rotate rotates the matrix by theta.
 func (g *GeoM) Rotate(theta float64) {
 	sin, cos := math.Sincos(theta)
-	g.Concat(GeoM{
-		elements: []float64{
-			cos, -sin, 0,
-			sin, cos, 0,
-		},
+	g.Concat(&GeoM{
+		a:      cos,
+		b:      -sin,
+		c:      sin,
+		d:      cos,
+		inited: true,
 	})
 }
 
diff --git a/vertices.go b/vertices.go
index 9eda99a8f..d66569067 100644
--- a/vertices.go
+++ b/vertices.go
@@ -31,13 +31,13 @@ func vertices(parts ImageParts, width, height int, geo *affine.GeoM) []float32 {
 	// TODO: This function should be in graphics package?
 	l := parts.Len()
 	vs := make([]float32, l*quadFloat32Num)
-	g := geo.UnsafeElements()
-	g0 := float32(g[0])
-	g1 := float32(g[1])
-	g2 := float32(g[3])
-	g3 := float32(g[4])
-	g4 := float32(g[2])
-	g5 := float32(g[5])
+	a, b, c, d, tx, ty := geo.Elements()
+	g0 := float32(a)
+	g1 := float32(b)
+	g2 := float32(c)
+	g3 := float32(d)
+	g4 := float32(tx)
+	g5 := float32(ty)
 	w := 1
 	h := 1
 	for w < width {