10-714：hw0

不过，我们要说的是，在Python版本中，应该大量使用numpy中的线性代数函数调用：尝试使用显式循环通常会使代码慢得多。

10-714课程的所有代码开发都可以在Google Colab环境中完成。但是，与其在colab笔记本内部广泛使用实际的代码块不同，你将大部分开发的代码编写为.py文件，这些文件将被自动下载到您的Google Drive中，你将主要使用笔记本运行shell脚本来测试和提交代码到自动评分系统（或者选择性地用于在开发过程中测试代码片段，但这不是必需的）。这是一种非常标准的Colab笔记本使用方法（通常人们更多地将其用作交互式编码环境，直接在笔记本中使用代码单元格）。然而，我们之所以按照这种方式使用它们的理由实际上很简单：除了是一个很好的基于云的笔记本环境之外，Colab还提供了非常好的访问“标准”基于云的GPU系统的方式，可以快速启动这些系统，从而使你能够开发一些后期（基于CUDA）的代码，而无需获得物理GPU的访问权限或自行设置CUDA库。话虽如此，你可以选择在任何您喜欢的环境中开发和提交代码，我们只是无法保证支持Colab以外的任何环境。

要开始，请通过从“文件”菜单中选择“在Drive中保存副本”来复制此笔记本文件，然后运行下面的代码块。这将将您的Google Drive文件夹加载到Colab

Q1

Q2

import gzip
import numpy as np
import struct

def parse_mnist(image_filename, label_filename):
    with gzip.open(image_filename, 'rb') as image_file, gzip.open(label_filename, 'rb') as label_file:
        # Read magic numbers from image and label files
        image_magic = struct.unpack('>I', image_file.read(4))[0]
        label_magic = struct.unpack('>I', label_file.read(4))[0]

        # Verify magic numbers
        if image_magic != 2051 or label_magic != 2049:
            raise ValueError("Invalid MNIST file format")

        # Read number of images and labels
        num_images = struct.unpack('>I', image_file.read(4))[0]
        num_labels = struct.unpack('>I', label_file.read(4))[0]

        # Verify that the number of images and labels match
        if num_images != num_labels:
            raise ValueError("Number of images and labels do not match")

        # Read image dimensions
        num_rows = struct.unpack('>I', image_file.read(4))[0]
        num_cols = struct.unpack('>I', image_file.read(4))[0]
        input_dim = num_rows * num_cols

        # Read images and labels into numpy arrays
        X = np.frombuffer(image_file.read(), dtype=np.uint8).reshape(num_images, input_dim)
        y = np.frombuffer(label_file.read(), dtype=np.uint8)

        # Normalize X to have values between 0.0 and 1.0
        X = X.astype(np.float32) / 255.0

        return X, y

• 数据集

在MNIST数据集中，训练集由两个文件组成：一个是包含标签的文件（train-labels-idx1-ubyte），另一个是包含图像数据的文件（train-images-idx3-ubyte）。

训练集标签文件（train-labels-idx1-ubyte）的格式描述如下：

• 偏移量0-3：32位整数，魔术数字（Magic Number），其值为0x00000801（2049），用于标识文件类型（高位优先）。

• 偏移量4-7：32位整数，60000，表示训练集中样本的数量。

• 偏移量8及之后：每个字节表示一个样本的标签，取值范围为0到9。

训练集图像文件（train-images-idx3-ubyte）的格式描述如下：

• 偏移量0-3：32位整数，魔术数字（Magic Number），其值为0x00000803（2051），用于标识文件类型。

• 偏移量4-7：32位整数，60000，表示训练集中图像的数量。

• 偏移量8-11：32位整数，28，表示图像的行数。

• 偏移量12-15：32位整数，28，表示图像的列数。

• 偏移量16及之后：每个字节表示一个像素值，取值范围为0到255。像素值以逐行的方式组织，每行的像素从左到右排列。

需要注意的是，图像的像素值为0到255，其中0表示背景（白色），255表示前景（黑色）

(**读取像素点的时候 ，是uint8 **)

• '>I'是一个格式化字符串，用于指定数据的解析格式。在这个上下文中，'>I'表示按照大端字节序（Big Endian）解析数据，并将其转换为无符号整数（unsigned integer）。

  • >表示使用大端字节序，即高位字节在前、低位字节在后的方式表示数据。

  • I表示将数据解析为一个无符号整数（32位）。

  通过使用struct.unpack('>I', data)，struct模块能够按照指定的格式解析二进制数据，并将其转换为对应的Python对象。在这种情况下，它将读取的4个字节数据解析为一个大端字节序的无符号整数。

• image_file.read()读取操作会从当前文件指针位置开始读取指定数量的字节，并且在读取后会更新文件指针的位置。每次调用read()方法会将文件指针向后移动读取的字节数。

• 对于图像数据，使用np.frombuffer(image_file.read(), dtype=np.uint8)将读取的文件内容转换为NumPy数组。frombuffer函数将字节流转换为数组，dtype=np.uint8指定了数组元素的数据类型为无符号8位整数，即像素值的范围为0-255。然后，使用.reshape(num_images, input_dim)将一维的数组变为二维数组

Q3

def softmax_loss(Z, y):
    """ Return softmax loss.  Note that for the purposes of this assignment,
    you don't need to worry about "nicely" scaling the numerical properties
    of the log-sum-exp computation, but can just compute this directly.

    Args:
        Z (np.ndarray[np.float32]): 2D numpy array of shape
            (batch_size, num_classes), containing the logit predictions for
            each class.
        y (np.ndarray[np.int8]): 1D numpy array of shape (batch_size, )
            containing the true label of each example.

    Returns:
        Average softmax loss over the sample.
    """
    ### BEGIN YOUR CODE
    batch_size=Z.shape[0]

    exp_Z=np.exp(Z)
    probs=exp_Z/np.sum(exp_Z,axis=1,keepdims=True)

    yprob=probs[np.arange(batch_size),y]



    loss=-np.log(yprob)

    ave_loss=np.mean(loss)

    return ave_loss

In this code, Z represents the logit predictions for each class, and y represents the true labels. Here's a breakdown of the steps:

1. Calculate the exponential of Z using np.exp(Z).

2. Compute the sum of the exponential values along the axis 1 (which represents the classes) using np.sum(exp_Z, axis=1, keepdims=True).

3. Divide each element of exp_Z by the sum to get the softmax probabilities for each example.

4. Select the probabilities for the correct class using advanced indexing probs[np.arange(batch_size), y], where np.arange(batch_size) generates an array [0, 1, ..., batch_size-1] to select the correct elements for each example.

5. Compute the negative logarithm of the correct probabilities using -np.log(correct_probs).

6. Finally, calculate the average loss over the batch by taking the mean of the log probabilities using np.mean(log_probs).

Q4

我好蠢啊，我才明白矩阵运算过来

$$abla_{\theta} l_{ce}(X\theta,y)=X^T(Z-I_y)$$

这样算出来的梯度是该批次的梯度和，所以需要除以$B$。所以矩阵运算其实就是做并行计算

跑一个epoch其实就是把数据集随机化以后，跑完全部的数据集

def softmax_regression_epoch(X, y, theta, lr = 0.1, batch=100):
    """ Run a single epoch of SGD for softmax regression on the data, using
    the step size lr and specified batch size.  This function should modify the
    theta matrix in place, and you should iterate through batches in X _without_
    randomizing the order.

    Args:
        X (np.ndarray[np.float32]): 2D input array of size
            (num_examples x input_dim).
        y (np.ndarray[np.uint8]): 1D class label array of size (num_examples,)
        theta (np.ndarrray[np.float32]): 2D array of softmax regression
            parameters, of shape (input_dim, num_classes)
        lr (float): step size (learning rate) for SGD
        batch (int): size of SGD minibatch

    Returns:
        None
    """
    ### BEGIN YOUR CODE
    num_example = X.shape[0]
    num_batch=num_example//batch
    num_class=theta.shape[1]




    for i in range(num_batch):
      X_batch=X[i*batch:(i+1)*batch]
      y_batch=y[i*batch:(i+1)*batch]
      Z=np.dot(X_batch,theta)
      probs=np.exp(Z)/ np.sum(np.exp(Z),axis=1,keepdims=True)
      I_y = np.zeros((batch, num_class))
      I_y[np.arange(batch), y_batch] = 1
      grad=np.dot(X_batch.T,(probs-I_y)) /batch
      theta-=lr*grad

Q5

def nn_epoch(X, y, W1, W2, lr = 0.1, batch=100):
    """ Run a single epoch of SGD for a two-layer neural network defined by the
    weights W1 and W2 (with no bias terms):
        logits = ReLU(X * W1) * W2
    The function should use the step size lr, and the specified batch size (and
    again, without randomizing the order of X).  It should modify the
    W1 and W2 matrices in place.

    Args:
        X (np.ndarray[np.float32]): 2D input array of size
            (num_examples x input_dim).
        y (np.ndarray[np.uint8]): 1D class label array of size (num_examples,)
        W1 (np.ndarray[np.float32]): 2D array of first layer weights, of shape
            (input_dim, hidden_dim)
        W2 (np.ndarray[np.float32]): 2D array of second layer weights, of shape
            (hidden_dim, num_classes)
        lr (float): step size (learning rate) for SGD
        batch (int): size of SGD minibatch

    Returns:
        None
    """
    ### BEGIN YOUR CODE
    num_example = X.shape[0]
    num_batch=num_example//batch
    num_class=W2.shape[1]

    for i in range(num_batch):
      X_batch=X[i*batch:(i+1)*batch]
      y_batch=y[i*batch:(i+1)*batch]
      Z1=relu(np.dot(X_batch,W1))
      Z1W2=np.dot(Z1,W2)
      probs=np.exp(Z1W2)/ np.sum(np.exp(Z1W2),axis=1,keepdims=True)
      I_y = np.zeros((batch, num_class))
      I_y[np.arange(batch), y_batch] = 1
      G2=probs-I_y
      G1=relu_derivative(Z1)*np.dot(G2,W2.T)
      grad1=np.dot(X_batch.T,G1)/batch
      grad2=np.dot(Z1.T,G2)/batch
      W1-=lr*grad1
      W2-=lr*grad2


    ### END YOUR CODE

Q6

$$Z=normalize(X\theta)$$

X=B乘n,$\theta$= n乘k,Z=B乘k

$$abla_{\theta} l_{ce}(X\theta,y)=\frac{1}{B} \bold{{X^T(Z-I_y)}}$$

每个batch

$$\theta=\theta-lr \cdot \nabla_{\theta} l_{ce}(X\theta,y)$$

参数

void softmax_regression_epoch_cpp(const float *X, const unsigned char *y,
								  float *theta, size_t m, size_t n, size_t k,
								  float lr, size_t batch)

void matrix_multi(const float *M,const float*N,float *P,size_t a,size_t b,size_t c){
    for (size_t row=0;row<a;row++){
        for (size_t col=0;col<c;col++){
            P[row*c+col]=0.0;
            for (size_t mid = 0; mid < b; mid++)
            {
                P[row*c+col]+=M[row*b+mid]*N[mid*c+col];
            } 
        }
    }
}


void normalize(float *M,size_t a, size_t b ){
    
    for (size_t row=0;row<a;row++){
        float max=0.0;
        float sum=0.0;
        for (size_t col=0;col<b;col++){
            max=std::max(max,M[row*b+col]);
        }

        for (size_t col=0;col<b;col++){
            M[row*b+col]=std::exp(M[row*b+col]-max);
            sum+= M[row*b+col];
        }

        for (size_t col=0;col<b;col++){
            M[row*b+col]/=sum;
        }
    }
}


void transpose(const float * M,float * N,size_t a,size_t b){
    for (size_t row=0;row<a;row++){
        for (size_t col=0;col<b;col++){
            N[col*a+row]=M[row*b+col];
        }
    }
}


void print(const float * M,size_t a,size_t b ){
    for (size_t row=0;row<a;row++){
        for (size_t col=0;col<b;col++){
            std::cout<<M[row*b+col]<<" ";
        }
        std::cout<<std::endl;
    }
    std::cout<<std::endl;
}




    void softmax_regression_epoch_cpp(const float *X, const unsigned char *y,
                                    float *theta, size_t m, size_t n, size_t k,
                                    float lr, size_t batch)
    {
        


        float * Z=new float[batch*k];//保存计算的logits
        float * grad=new  float[n*k];//保存计算的梯度
        float * X_T=new float[n*batch];
      
      
      	

        for (size_t  batch_start=0;batch_start<m;batch_start+=batch)
        {
            size_t batch_end=std::min(m,batch_start+batch);
          
          
          	//Z=normalize(X\theta)
            matrix_multi(X+batch_start*n,theta,Z,batch_end-batch_start,n,k);
          
          
            normalize(Z,batch_end-batch_start,k);
         
          
          	//Z-I_y
            for (size_t i = 0; i <(batch_end - batch_start); i++)
            {
                    Z[i*k+y[i + batch_start]] -= 1.0;
            }
          
  					//X_T
            transpose(X+batch_start*n,X_T,batch_end-batch_start,n);
          
          
          	//grad=X_T *(Z-I_y)
            matrix_multi(X_T,Z,grad,n,batch_end-batch_start,k);

            for (size_t i = 0; i < n*k; i++)
            {
                theta[i]-=lr*grad[i]/(batch_end-batch_start);
            }
            
        }

        delete[] Z;        
        delete[] grad;
        delete[] X_T;

       
    }

why

softmax函数的目标是将一个k维的向量映射到另一个k维的向量，其值在0和1之间且总和为1，可以被解释为概率。softmax函数定义如下：

$$softmax(x_i) = \frac{e^{x_i}}{\sum e^{x_j}}$$

其中，i和j都在1到k的范围内。

但是，直接这样计算可能会有数值稳定性问题。当x的某个值非常大的时候，${e^{x_i}}$可能会超出浮点数的表示范围，导致数值溢出。

为了解决这个问题，一种常见的技巧是在计算softmax的时候首先减去x中的最大值。即：

$$softmax(x_i) = \frac{e^{x_i-max(x)}}{\sum e^{x_j-max(x)}}$$

可以证明这样做并不改变softmax的结果，但是可以防止因为x的值过大而导致的数值溢出。