diff --git a/python/binding.cpp b/python/binding.cpp
index 82ad36d5..e3e18224 100644
--- a/python/binding.cpp
+++ b/python/binding.cpp
@@ -486,7 +486,7 @@ NB_MODULE(scaluq_core, m) {
     DEF_TWO_QUBIT_GATE(TwoQubitMatrixGate, "Specific class of double-qubit dense matrix gate.")
         .def(
             "matrix",
-            [](const TwoQubitMatrixGate &gate) { gate->matrix(); },
+            [](const TwoQubitMatrixGate &gate) { return gate->matrix(); },
             "Get property `matrix`.");
 
     DEF_GATE(FusedSwapGate,
@@ -555,6 +555,7 @@ NB_MODULE(scaluq_core, m) {
     DEF_GATE_FACTORY(U1);
     DEF_GATE_FACTORY(U2);
     DEF_GATE_FACTORY(U3);
+    DEF_GATE_FACTORY(OneQubitMatrix);
     DEF_GATE_FACTORY(CX);
     mgate.def("CNot",
               &gate::CX,
@@ -562,8 +563,14 @@ NB_MODULE(scaluq_core, m) {
     DEF_GATE_FACTORY(CZ);
     DEF_GATE_FACTORY(Swap);
     DEF_GATE_FACTORY(FusedSwap);
+    DEF_GATE_FACTORY(TwoQubitMatrix);
     DEF_GATE_FACTORY(Pauli);
     DEF_GATE_FACTORY(PauliRotation);
+    mgate.def("DenseMatrix",
+              &gate::DenseMatrix,
+              "Generate general Gate class instance of DenseMatrix. IGate, OneQubitMatrixGate or "
+              "TwoQubitMatrixGate correspond to len(target) is created. The case len(target) >= 3 "
+              "is currently not supported.");
     DEF_GATE_FACTORY(Probablistic);
 
     nb::enum_<ParamGateType>(m, "ParamGateType", "Enum of ParamGate Type.")
diff --git a/scaluq/gate/gate_factory.hpp b/scaluq/gate/gate_factory.hpp
index f305575c..3728ead7 100644
--- a/scaluq/gate/gate_factory.hpp
+++ b/scaluq/gate/gate_factory.hpp
@@ -84,6 +84,9 @@ inline Gate U2(UINT target, double phi, double lambda) {
 inline Gate U3(UINT target, double theta, double phi, double lambda) {
     return internal::GateFactory::create_gate<internal::U3GateImpl>(target, theta, phi, lambda);
 }
+inline Gate OneQubitMatrix(UINT target, const std::array<std::array<Complex, 2>, 2>& matrix) {
+    return internal::GateFactory::create_gate<internal::OneQubitMatrixGateImpl>(target, matrix);
+}
 inline Gate CX(UINT control, UINT target) {
     return internal::GateFactory::create_gate<internal::CXGateImpl>(control, target);
 }
@@ -98,12 +101,58 @@ inline Gate FusedSwap(UINT qubit_index1, UINT qubit_index2, UINT block_size) {
     return internal::GateFactory::create_gate<internal::FusedSwapGateImpl>(
         qubit_index1, qubit_index2, block_size);
 }
+inline Gate TwoQubitMatrix(UINT target1,
+                           UINT target2,
+                           const std::array<std::array<Complex, 4>, 4>& matrix) {
+    return internal::GateFactory::create_gate<internal::TwoQubitMatrixGateImpl>(
+        target1, target2, matrix);
+}
 inline Gate Pauli(const PauliOperator& pauli) {
     return internal::GateFactory::create_gate<internal::PauliGateImpl>(pauli);
 }
 inline Gate PauliRotation(const PauliOperator& pauli, double angle) {
     return internal::GateFactory::create_gate<internal::PauliRotationGateImpl>(pauli, angle);
 }
+inline Gate DenseMatrix(const std::vector<UINT>& targets, const ComplexMatrix& matrix) {
+    UINT nqubits = targets.size();
+    UINT dim = 1ULL << nqubits;
+    if (static_cast<UINT>(matrix.rows()) != dim || static_cast<UINT>(matrix.cols()) != dim) {
+        throw std::runtime_error(
+            "gate::DenseMatrix(const std::vector<UINT>&, const ComplexMatrix&): matrix size must "
+            "be 2^{n_qubits} x 2^{n_qubits}.");
+    }
+    if (targets.size() == 0) return I();
+    if (targets.size() == 1) {
+        return OneQubitMatrix(targets[0],
+                              std::array{std::array{Complex(matrix(0, 0)), Complex(matrix(0, 1))},
+                                         std::array{Complex(matrix(1, 0)), Complex(matrix(1, 1))}});
+    }
+    if (targets.size() == 2) {
+        return TwoQubitMatrix(targets[0],
+                              targets[1],
+                              std::array{
+                                  std::array{Complex(matrix(0, 0)),
+                                             Complex(matrix(0, 1)),
+                                             Complex(matrix(0, 2)),
+                                             Complex(matrix(0, 3))},
+                                  std::array{Complex(matrix(1, 0)),
+                                             Complex(matrix(1, 1)),
+                                             Complex(matrix(1, 2)),
+                                             Complex(matrix(1, 3))},
+                                  std::array{Complex(matrix(2, 0)),
+                                             Complex(matrix(2, 1)),
+                                             Complex(matrix(2, 2)),
+                                             Complex(matrix(2, 3))},
+                                  std::array{Complex(matrix(3, 0)),
+                                             Complex(matrix(3, 1)),
+                                             Complex(matrix(3, 2)),
+                                             Complex(matrix(3, 3))},
+                              });
+    }
+    throw std::runtime_error(
+        "gate::DenseMatrix(const std::vector<UINT>&, const ComplexMatrix&): DenseMatrix gate more "
+        "than two qubits is not implemented yet.");
+}
 inline Gate Probablistic(const std::vector<double>& distribution,
                          const std::vector<Gate>& gate_list) {
     return internal::GateFactory::create_gate<internal::ProbablisticGateImpl>(distribution,
diff --git a/scaluq/gate/gate_matrix.hpp b/scaluq/gate/gate_matrix.hpp
index a02a68b1..53b6500c 100644
--- a/scaluq/gate/gate_matrix.hpp
+++ b/scaluq/gate/gate_matrix.hpp
@@ -57,8 +57,8 @@ class TwoQubitMatrixGateImpl : public TwoQubitGateBase {
                            UINT target2,
                            const std::array<std::array<Complex, 4>, 4>& matrix)
         : TwoQubitGateBase(target1, target2) {
-        for (UINT i : std::views::iota(4)) {
-            for (UINT j : std::views::iota(4)) {
+        for (UINT i : std::views::iota(0, 4)) {
+            for (UINT j : std::views::iota(0, 4)) {
                 _matrix.val[i][j] = matrix[i][j];
             }
         }
@@ -66,8 +66,8 @@ class TwoQubitMatrixGateImpl : public TwoQubitGateBase {
 
     std::array<std::array<Complex, 4>, 4> matrix() const {
         std::array<std::array<Complex, 4>, 4> matrix;
-        for (UINT i : std::views::iota(4)) {
-            for (UINT j : std::views::iota(4)) {
+        for (UINT i : std::views::iota(0, 4)) {
+            for (UINT j : std::views::iota(0, 4)) {
                 matrix[i][j] = _matrix.val[i][j];
             }
         }
@@ -77,8 +77,8 @@ class TwoQubitMatrixGateImpl : public TwoQubitGateBase {
     Gate copy() const override { return std::make_shared<TwoQubitMatrixGateImpl>(*this); }
     Gate get_inverse() const override {
         std::array<std::array<Complex, 4>, 4> matrix_dag;
-        for (UINT i : std::views::iota(4)) {
-            for (UINT j : std::views::iota(4)) {
+        for (UINT i : std::views::iota(0, 4)) {
+            for (UINT j : std::views::iota(0, 4)) {
                 matrix_dag[i][j] = Kokkos::conj(_matrix.val[j][i]);
             }
         }